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simplification  is  discussed  in  the  thesis  of  CAPT  Kevin  D. 
Jenkins,  U.S.  Marine  Corps.  In  this  thesis  we  first  discuss 
a  computational  investigation  of  two  methods  by  which  homotopy 
classes  can  be  named.  Next,  a  computational  heuristic  is 
presented  that  finds  the  lower  bound  for  a  path  in  a  class. 
Finally,  the  true  shortest  path  is  found  by  searching  these 
classes  in  order  of  increasing  lower  bound.  One  application 
of  this  study  is  in  the  area  of  robotic  path  planning. 


111 


L.I 


THESIS  DISCLAIMER 

The  reader  is  cautioned  that  computer  programs  developed 
in  this  research  may  not  have  been  exercised  for  all  cases  of 
interest.  While  every  effort  has  been  made,  within  the  time 
available,  to  ensure  that  the  programs  are  free  of 
computational  and  logic  errors,  they  cannot  be  considered 
validated.  Any  application  of  these  programs  without 
additional  verification  is  at  the  risk  of  the  user. 
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I .     INTRODUCTION 

The  industrial  community  relies  today,  in  increasing 
measure,  upon  robots  to  handle  a  multitude  of  tasks.  A  need 
remains  to  develop  robots,  and  fixed  robot  arms,  with  the 
ability  to  roam  freely  among  obstacles.  In  this  regard,  the 
calculation  of  shortest  paths  is  of  obvious  importance. 

A.  THE  PROBLEM 

As  robot  arms  move,  they  must  vary  their  configurations  to 
position  their  end-effector  or  "hand".  Once  obstacles  are 
introduced,  freedom  of  movement  of  the  manipulator  may  become 
drastically  reduced.  The  problem  considered  here  is  one  of 
finding  the  shortest  path  from  a  starting  point  to  a 
destination  point  in  a  planar  region  containing  obstacles. 
This  paper  addresses  several  parts  of  the  solution  to  this 
problem. 

B.  OVERVIEW  OF  THE  SOLUTION 

Between  any  two  points  in  the  plane,  a  path  joining  them 
may  be  chosen  from  an  uncountably  infinite  set  of 
alternatives.  It  is  desirable  to  choose  the  shortest  path 
from  this  set  which,  due  to  the  presence  of  obstacles,  may  not 
be  a  straight  line.  An  exhaustive  search  of  the  collection  of 
possible  paths  is  impossible,  so  another  method  must  be 
developed. 


The  set  of  possible  paths  is  partitioned  into  a  collection 
of  equivalence  classes—mutual  ly  exclusive  subsets  which 
collectively  exhaust  the  partitioned  set.  This  partition  is 
induced  when  an  equivalence  relation  is  defined  on  the  set  of 
possible  paths.  This  set  of  equivalence  classes  is  countably 
infinite. 

Naming  conventions  are  established  which  associate  with 
each  class  a  character  string  which  allows  the  classes  to  be 
represented  in  a  computer. 

Next,  a  finite  list  of  candidate  classes  is  produced  to 
search. 

A  heuristic  is  then  applied  to  this  finite  list  of  classes 
to  place  them  in  an  order  which  saves  computational  effort. 

The  final  step  begins  with  this  ordered  list  of  candidate 
classes  and  solves  the  shortest  path  problem  class  by  class  in 
listed  order.  Savings  of  computational  effort  are  realized 
when  the  class-by-class  solution  process  can  terminate  without 
exhausting  the  ordered  list  of  candidate  classes. 

The  key  ideas  in  the  above  overview  are  now  considered  in 
somewhat  greater  detail. 

C.   THE  APPROACH 

1.   The  Topology 

Let  P  be  the  uncountably  infinite  set  of  all 
continuous,  obstacle-avoiding  paths  from  a  starting  point  to 
a  destination  point,  denoted  a  and  z  respectively.   If  p^  and 


pj  are  two  paths  in  P,  which  have  the  same  starting  point  and 
destination  point,  we  say  that  p:  is  homotopic  to  p^  if  p-  can 
be  mapped  to  pj  under  a  continuous  function  (with  both 
endpoints  fixed  in  place)  without  encroaching  on  any  obstacles 
[Ref.  l:p.  223].  Clearly,  the  homotopy  relation  is  reflexive, 
symmetric  and  transitive  and  therefore  defines  an  equivalence 
relation  [Ref.  l:p.  223].  This  relation  induces  a  partition 
of  P  into  a  countably  infinite  collection  of  equivalence 
classes,  known  as  homotopy  classes  [Ref.  l:p.  223]. 
2.   Establishing  the  Naming  Convention 

In  order  to  name  homotopy  classes,  a  reference  frame 
is  established  to  represent  the  topological  relationships  in 
the  plane  with  obstacles.  For  a  given  path,  p,  a  string  of 
characters,  R(p),  is  recorded  which  encodes  information 
concerning  the  relationship  of  the  path  to  the  obstacles. 

Two  algorithms  are  then  presented  which  accept  R(p)  as 
input.  These  algorithms  have  the  following  important  property 
which  implicitly  defines  names  for  the  homotopy  classes:  If 
p^  and  p.  are  coterminal  paths  and  R(pp  and  R(p^)  are  input  to 
either  of  the  two  algorithms,  then  the  outputs  are  identical 
if  and  only  if  p^  is  homotopic  to  p^  [Ref.  2]. 

Chapter  III  presents  a  computational  investigation 
which  closely  examines  the  methods  employed  by  these  two 
algorithms  to  produce  the  names  for  these  homotopy  classes. 


3.  Generation  of  Names  for  Candidate  Equivalence  Classes 

An  edge-labeled  graph  is  constructed  which  models  the 
topological  relationships  within  the  region,  where  the  nodes 
represent  subregions  induced  by  the  reference  frame.  Nodes 
are  connected  if  their  corresponding  subregions  are  adjacent. 
The  names  of  the  desired  homotopy  classes  may  be  obtained  by 
traversing  this  graph.  This  traversal  produces  a  list  of 
candidate  classes  which  contain  all  paths  of  minimal  length. 

The  location  of  the  origin  affects  the  manner  in  which 
the  plane  is  divided  into  wedges.  In  this  regard,  the  graph 
which  is  created,  given  the  location  of  the  origin  in  one 
instance,  may  not  be  the  same  as  the  graph  obtained  if  the 
origin  is  moved  to  a  new  point. 

4.  Heuristic  Ordering  of  Candidate  Homotopy  Classes 
For  each  class  on  the  list,  a  lower  bound  on  the 

length  of  its  shortest  representative  path  is  constructed. 
The  list  of  classes  is  then  arranged  into  increasing  order  of 
these  bounds.  To  obtain  the  bounds,  a  point  is  first  chosen 
within  each  obstacle.  A  contraction  deformation  is  then  be 
applied  to  each  obstacle  to  "shrink"  the  obstacle  to  that 
chosen  point. 

Now  some  class  is  fixed  and  its  shortest  path  is 
examined.  As  all  obstacles  are  simultaneously  contracted  to 
their  representing  points,  this  shortest  path  has  a  limiting 
position  which  is  polygonal.  The  length  of  this  polygonal 
path  is  the  lower  bound  which  is  associated  with  that  class. 


This  length  can  be  calculated  from  the  class  name  and  the 
representing  points  without  explicitly  defining  any 
contraction  deformation.  This  thesis  presents  this  method  to 
determine  a  bound  for  the  shortest  path  in  a  given  class. 

5.  Class  by  Class  Solution  of  the  Shortest  Path  Problem 
In  the  final  step  of  the  solution,  the  classes  on  the 
ordered  list  of  candidates  are  considered.  The  first  class  is 
removed  from  the  list--that  class  with  the  smallest  associated 
bound.  The  true  shortest  path  in  the  named  class  is  found  by 
reversing  the  contraction  previously  applied  to  the  obstacles, 
thereby  transforming  the  polygonal  path  whose  length  provided 
the  lower  bound  into  the  true  shortest  path.  If  the  length  of 
this  path  is  smaller  than  the  bound  associated  with  the  class 
on  the  top  of  the  remaining  list,  the  search  is  stopped. 
Otherwise,  the  first  class  is  removed  from  that  remaining  list 
and  the  procedure  is  repeated.  This  procedure,  which  is 
contained  in  the  thesis,  continues  until  the  condition 
specified  above  is  met,  and  the  shortest  path  has  been  found. 

D.   Summary 

The  solution  to  the  problem  of  searching  for  the  shortest 
path  between  two  points  in  the  plane  with  obstacles  begins 
with  consideration  of  a  set  of  paths  which  is  uncountably 
infinite.  Through  the  homotopy  relation,  this  set  is  reduced 
to  a  countably  infinite  set.  The  new  set  is  further  reduced 
to  a  finite  list  by  modeling  the  region  with  a  graph  and 


applying  an  optimality  principle.  This  final  list  contains 
only  those  homotopy  classes  containing  paths  which  are  not 
sel f -crossing .   From  this  list  a  shortest  path  is  found. 

Computational  effort  is  further  reduced  through  the  use  of 
a  heuristic  which  orders  the  list  of  candidate  classes  by 
increasing  order  of  their  lower  bounds.  The  heuristic  used 
also  facilitates  a  methodical  search  of  the  classes  while 
solving  the  shortest  path  problem.  The  use  of  this  heuristic 
does  not,  however,  imply  that  the  solution  is  approximate. 
The  solution  to  the  problem  will  be  exact  using  this  method. 

E.   THE  CONTRIBUTION 

Chapter  II  of  this  thesis  presents  the  methods  used  to 
establish  a  reference  frame  given  the  plane  with  obstacles  and 
an  algorithm  which  is  used  to  generate  the  homotopy  class  name 
of  a  given  path. 

Chapter  III  presents  a  computational  investigation  which 
addresses  the  question  of  whether  or  not  two  algorithms  which 
are  used  to  name  homotopy  classes  actually  provide  the  same 
results . 

Chapter  IV  introduces  a  method  to  find  a  lower  bound  for 
the  shortest  path  in  a  class. 

Chapter  V  shows  how  to  find  the  true  shortest  path  in  a 
class . 


II.   NAMING  EQUIVALENCE  CLASSES 

A.   INTRODUCTION 

In  this  chapter,  a  labelling  scheme  is  established  to 
represent  the  equivalence  classes.  This  notation  is  used 
throughout  the  paper  to  organize  the  computational  search  for 
the  shortest  path. 

An  algorithm  is  presented  that  for  a  character 
representation  of  a  path,  p,  names  the  path's  respective 
equivalence  class.  It  can  be  shown  that,  after  being 
processed  by  the  algorithm,  different  character 
representations  of  paths  p^  and  p^  yield  the  same  output 
exactly  when  and  only  when  the  two  paths  are  in  the  same 
equivalence  class  [Ref.  2].  Thus,  the  name  of  a  class  will  be 
the  string  obtained  by  applying  the  algorithm  to  any 
representative  in  the  class. 

The  computation  of  class  names  depends  on  a  reference 
frame  which  in  turn  depends  on  a  collection  of  obstacles. 
Although  more  than  one  reference  frame  can  be  drawn  for  a 
particular  collection  of  obstacles,  the  choice  of  a  particular 
frame  fixes  the  representation  of  all  homotopy  classes. 

Once  a  reference  frame  is  developed,  the  name  of  the 
homotopy  class  for  a  path  can  be  determined  by  a  two-step 
procedure.  First,  a  character  string  R(p)  is  calculated  which 
encodes  certain  information  about  the  path  taken  through  the 


obstacles.  Second,  this  character  string  is  accepted  as  input 
to  an  algorithm  which  then  produces  a  name  for  the  equivalence 
class  to  which  the  particular  path  belongs.  This  algorithm 
produces  the  class  names  in  terms  of  the  same  alphabet  used  to 
create  the  initial  character  string.  The  results  from  the 
algorithm  described,  Algorithm  1,  are  used  throughout  this 
analysis  as  the  class  name  associated  with  a  given  path. 

B.   CONSTRUCTION  OF  REFERENCE  FRAME 

Let  b^  be  an  arbitrary  point  in  obstacle  B^,  k  =  1,2,  .  .  .  ,n, 
where  n  is  the  number  of  obstacles  in  the  region.  A  point  c 
is  chosen  and  a  line  drawn  through  each  b^,  infinite  in  extent 
in  each  direction  and  having  the  following  properties:  there 
is  an  open  disk,  6,  centered  at  c  such  that  506^=0  for  all 
k=l,...,n,  and  the  n  lines  connecting  c  with  the  points  b^  are 
distinct.  Such  a  point  c  can  always  be  found  as  the  above  two 
conditions  are  satisfied  by  any  point  in  the  planar  region 
that  is  neither  on  an  obstacle  nor  on  the  n(n-l)/2  lines 
determined  by  pairwise  choices  of  distinct  b^. 

To  draw  a  reference  frame,  the  n  lines  are  first 
constructed  joining  c  with  each  b^ .  The  line  from  c  to  b^  is 
labeled  as  L^ .  Each  line  is  then  partitioned  into  two 
directed,  semi-infinite  rays  and  one  finite  length  line 
segment.  The  ray  directed  from  c  in  the  direction  away  from 
b^  is  called  o^ .  The  ray  emanating  from  b^  and  away  from  c  is 
denoted  B^.    The  remaining  line  segment,   [c^b^],  is  also 
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denoted  at..  The  reference  frame  is  the  collection  of  line 
segments  and  rays  so  constructed,  as  illustrated  in  Figure 
2.1. 


Figure  2.1   A  Reference  Frame  With  n=2  Obstacles 


C.   RAW  CHARACTER  STRINGS 

A  reference  frame  is  fixed  in  a  region  T.  Let  a  and  z  be 
points  in  T  and  let  p  be  a  directed  path  in  T  from  a  to  z  . 
Then  the  raw  string  of  p,  denoted  R(p),  is  defined  to  be  the 
ordered  sequence  of  characters  obtained  by  following  p  from  a 
to  z  and  recording,  in  order,  the  names  of  the  rays  that  are 
crossed . 

Two  special  cases  must  be  addressed  so  as  to  make  the 
above  definition  complete.  First,  in  the  case  that  p  crosses 
no  reference  rays,  we  let  R(p)  =  e,  where  e  denotes  the  empty 
string.    Second,   if  p  crosses   through  c   (simultaneously 


crossing  all  a.)  the  names  of  all  a:  will  be  recorded  in  order 

j  j 

of  increasing  subscript . 

Thus,  raw  strings  are  of  the  form 

R(p)  =  Xlx2.  .  .xB 


where  each  X;  belongs  to  the  alphabet 
A  =  {e,  o^,  Q2  '  •  •  •  '  an'  Cl 


'  Bn> 


;th 


The  character  X:  ( j  =  1,  2,  .  .  .  ,  m)  is  the  name  of  the  j 
reference  ray  crossed  by  p.  Figure  2.2  shows  a  pair  of  paths 
connecting  the  two  points  a  and  z  in  the  region  with  two 
obstacl es . 


Figure  2.2   Paths  a,Bi  and  a,a 


2"1 


:lu2 


D.   ALGORITHM  1 

Algorithm  1  accepts  as  input  a  raw  string  R(p)  for  any 
path  p  in  T  and  any  reference  frame.  The  output  is  a  string 
denoted  C(R(p))  of  characters  also  chosen  from  the  alphabet  A. 
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This  output  string  C(R(p))  is  the  name  of  the  homotopy  class 
to  which  p  belongs. 

The  algorithm  is  presented  in  terms  of  two  functions.  The 
first  is  the  sorting  function  o.  Let  S  =  Xj...xffl  be  any  string 
over  A.  If  S  contains  a  two  character  substring  XiXim  =  ct'Ct' 
with  i<j,  then  o(S)  is  the  string  which  results  by  reversing 
the  order  of  the  leftmost  such  substring.  Figure  2.3  depicts 
two  such  strings  and  makes  clear  that  such  paths  are 
homotopic.  If  S  contains  no  such  two  character  substring, 
then  o(S)  =  S. 


Figure  2.3   Homotopic  Alpha  Strings  <^2a[    an^  a\a2 


Repeated  application  of   o  sorts  all   substrings   that 

consist   entirely   of   ctj's   into   non-decreasing   order   of 

j 

subscript.   2(S)  is  defined  to  be  the  string  which  results 
when  this  sorting  is  complete. 

The  second  function  is  the  cancellation  function  %.        If 
string  S  contains  a  two  character  substring  x^x^  with  x^  =  Xj^j 
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then  x(s)  i-s  the  string  which  results  by  removing  the  leftmost 
such  two  character  substring.  Otherwise  x(s)  =  s-  Such 
cancellation  is  intuitive,  for  a  reference  ray  that  is  crossed 
twice  consecutively  is  equivalent  to  one  that  is  not  crossed 
at  all.  Thus  repeated  application  of  %  cancels  all  pairs  of 
adjacent  like  occurrences  of  identical  characters.  Let  X(S) 
be  defined  as  the  string  which  results  when  all  such  possible 
cancellation  is  complete. 

With  these  definitions  complete,  Algorithm  1  is  given  in 
Figure  2.4. 


begin 
k  -  0 
Sk  -  R(p) 

WHILE  Sk  not  equal  to  X(Z(Sk)) 
k  -  k+1 


sk  -  x(z(sk.1)) 

WHT" 


END  WHILE 
C(R(p))  -  Sk 
END 


Figure  2.4   Algorithm  1 

The  output  from  Algorithm  1,  C(R(p)),  is  called  the 
canonical  representation  of  p  and  is  the  unique  name  for  the 
homotopy  class  of  p  [Ref.  2].  That  name  is  used  throughout 
this  paper. 
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III.   A  COMPUTATIONAL  INVESTIGATION 

A.  INTRODUCTION 

This  chapter  describes  a  computational  investigation  which 
was  conducted  as  an  aid  to  verifying  that  the  homotopy  class 
names  produced  by  Algorithm  1  are  equivalent  to  those  produced 
by  a  known  method  of  naming  homotopy  classes  [Ref.  2]. 

A  second  algorithm,  Algorithm  2,  is  introduced  here  to 
forge  a  link  between  the  class  names  of  Chapter  II  and  the 
well-known  fundamental  group.  Two  processes  are  presented  to 
exploit  this  link:  If  Algorithm  1  truly  names  the  classes, 
then  both  procedures  should  produce  the  same  output  for  every 
path  tested.  One  million  test  cases  were  examined  and  no 
counter-examples  were  found.  While  this  computational 
evidence  is  not  a  proof,  it  did  provide  initial  support  before 
the  proof  of  the  claim  that  Algorithm  1  names  homotopy  classes 
was  found. [Ref.  2] 

B.  ALGORITHM  2 

1.   Fundamental  Group 

The  class  name  obtained  from  Algorithm  1  was  written 
in  terms  of  the  alphabet  A  defined  in  Chapter  II.  We  present 
an  algorithm  here  where  the  name  obtained  is  expressed  in 
terms  of  the  fundamental  group  of  a  topological  space  T. 
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The  basic  idea  for  producing  a  fundamental  group  for 
this  problem  is  to  regard  the  paths  in  T  as  elements  of  the 
group  and  path  concatenation  (*)  as  the  group  operation. 
There  are  two  problems  with  considering  the  paths  as  the  group 
elements.  First,  it  is  not  necessarily  possible  to 
concatenate  any  two  paths  in  T.  In  order  to  concatenate  two 
paths  using  the  operation  (*),  the  first  path  must  end  at  the 
point  where  the  second  begins.  In  this  analysis  the  test 
paths  all  begin  and  end  at  the  point  a.  This  point  is  called 
a  base  point  and  concatenation  is  possible  for  every  pair  of 
paths . 

A  second  issue  is  that  an  inverse  is  not  well  defined 
for  paths.  If  the  identity  is  defined  as  the  empty  string, 
i.e.,  stay  at  the  base  point,  then  when  a  path  and  its  inverse 
are  concatenated  we  have  traced  out  a  path  which  is  not  equal 
to  the  identity.  To  avoid  this  problem,  the  homotopy  classes 
of  paths  are  considered.  Then  all  paths  that  have  an  empty 
character  string  representation  will  be  in  the  identity  class. 
Also,  any  class  followed  by  its  inverse  will  equal  the 
identity.  In  Figure  3.1,  g^  represents  the  class  ol^Q^,  and 
(g^)  represents  the  class  B^a^ .  So,  g^  *  (gj)  represents  the 
class  ctjfijBjaj  which  reduces  to  the  identity  class  by  the  rules 
of  Algorithm  1. 

The  group  elements,  therefore,  are  homotopy  classes  of 
loops  around  obstacles  based  at  a  common  point  and  the  group 
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operation  is  defined  in  terms  of  concatenation  (*)  of  these 
cl asses . 


g   :  class   a  B 

g   :  class   p1a1 

g  *  g~  :  PjCCjOjPp  e 


Figure  3.1   A  Homotopy  Class  and  its  Inverse 


Given  a  base  point  a,  we  let  p  be  any  loop  beginning 
and  ending  at  a.  We  let  [p^]  denote  the  homotopy  class  of  p^  . 
So  the  set  of  group  elements  is  {[p^]  such  that  p^  is  a  member 
of  T,  where  p^  is  a  loop  based  at  a}.  This  set  will  be  called 
G  with  the  elements  denoted  g>.  Figure  3.2  illustrates 
several  paths  and  their  homotopy  classes  given  n=2  obstacles. 


3 


Figure    3.2       [pj]    =    [p2]    *    [p3] 
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Then  the  group  operation  (*)  can  be  defined  by 

[p]  *  [q]  =  [pq]. 
It  is  important  to  note  that  the  fundamental  group  is 
finitely  generated.  As  Figure  3.3  illustrates,  the  generators 
are  not  unique.   However,  by  fixing  a  set  of  generators,  the 
class  names  become  fixed. 


2 


2 


Figure  3.3   Alternate  Generators  For  Fundamental  Group  of  the 

Space  With  Two  Obstacles 


With  the  above  information  describing  the  fundamental 
group  representation  and  a  few  necessary  definitions, 
Algorithm  2  is  given  below. 
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2.   Functions  in  Algorithm 

a.  Side  Array 

Consider  the  reference  lines  L^  constructed  in 
Chapter  II  to  be  oriented  rays  with  direction  from  c  to  b^.  It 
then  becomes  reasonable  to  discuss  a  'right'  and  a  'left'  side 
of  those  lines.  A  moving  point  is  allowed  to  trace  a  path  p 
from  start  point  to  destination  point.  Side(k)  is  a  function 
which  defines  the  side  of  L^  on  which  the  moving  point  lies. 
The  output  is  either  'left'  or  'right'.  The  output  is  never 
'on'  because  this  routine  is  used  only  after  a  complete 
crossing  of  L^  takes  place.  A  crossing  is  considered  to  be 
complete  when  the  moving  point  leaves  L^  to  one  side  after 
having  met  L^  from  the  opposite  side. 

b.  Switch  Function 

The  switch  function  is  defined  by: 

(    'right'  if  side  (k)  =  'left' 
switch  (side(k))  =  { 

I    'left'  if  side  (k)  =  'right' 

As  the  moving  point  traces  p,  each  time  a  reference  line  L^  is 

crossed,  the  switch  function  will  be  applied  to  indicate  on 

which  side  of  L^  the  moving  point  lies. 

c.  Index  Function 

Let  x^  be  any  character  from  the  alphabet  A 
representing  the  kth  element  in  R(p).  Let  indexCx^)  =  j  if  Xi 
equals  ctj  or  Ry  Thus,  the  k  crossing  completed  by  the  moving 
point  is  a  crossing  of  Lj  .   So  if  the  moving  point  is  on  the 
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left  side  of  Lj  initially  and  this  line  is  now  crossed  on 
either  side  of  bj,  then  side(3)  would  equal  'right'. 
d.   The  Algorithm 

Using  the  functions  described  above,  Algorithm  2 
is  given  in  Figure  3.4.  The  algorithm  works  in  two  phases. 
The  first  phase  initializes  the  array  side(k)  for  k  =  1, . . . ,n. 
The  second  phase  reads  the  raw  string  from  left  to  right  and 
adds  an  element  to  the  fundamental  group  representation  for 
each  fi  crossing.  When  the  end  of  the  input  string  is  reached, 
the  output  is  a  shorter  string  with  one  character 
corresponding  to  every  13  element  in  the  original  raw  string. 
Each  character  represents  a  generator  or  its  inverse. 


begin 

input  R(p)  =  XjXj  .  .  .xm 

for  j  =  1 ,  .  .  .  ,  n 

If  a  is  to  right  of  reference  line  L.  then 

side  (j)  -  right 
else 

side  (j)  «-  left 
end  if 
end  for 
G  -  [a] 

for  k  =  1 , . . . ,  m 
i  =  index(x^) 

side(i)  ♦-  switch(side(i  ) ) 
if  Xl  =  C„  for  some  r,  then 

if  sideCindexCx^) )  =  left,  then 

G  -  G  *  gi 
else 

G   -  G   *    (g[)'1 
end  if 
end  if 
end  for 
F(R(p))  -  G 
end 

Figure  3.4   Algorithm  2 
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e.   Fundamental  Group  Cancellation  Function 

The  cancellation  rules  in  the  fundamental  group 
differ  somewhat  from  those  for  the  raw  string.  Although  sets 
of  generators  are  not  unique,  we  can  obtain  a  unique 
representation  of  each  class  with  respect  to  a  particular  set 
of  generators.  For  any  given  set  of  generators,  every 
homotopy  class  can  be  represented  as  a  product  of  these 
generators  and  their  inverses.  Even  so,  this  representation 
is  only  unique  after  cancellation  is  applied.  The 
cancellation  rule  follows. 

Let  G  =  {g^,  .  .  .  ,  gn)  be  a  set  of  generators  of  the 
fundamental  group  of  T  (base  point  a)  and  let  Y  =  y^.-.yj,  be 
a  representation  of  some  homotopy  class  in  terms  of  the  g^  and 
their  inverses,  i.e.  for  each  i  =  1,...,  m,  y^  =  g:  or  (g^) 
for  some  j,k  =  1,...,  n. 

The  cancellation  function  k  is  defined  as  follows. 
If  Y  contains  a  two  character  substring  y^y^i  with  y^  =  (y^j) 
then  k(Y)  is  the  string  which  results  by  removing  both  y^  and 
y^j,  in  their  leftmost  occurrence,  from  Y.  Finitely  many 
repeated  applications  of  k  produce  a  string  in  which  no 
further  cancellation  is  possible.  We  define  this  string  as 
K(Y). 
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C.   THE  COMPUTATIONAL  INVESTIGATION 
1.   The  Approach 

A  graphical  representation  of  the  two  procedures  is 
provided  in  Figure  3.5.  The  investigation  accepts  a  random 
path  p  as  input  and  generates  two  names  for  the  homotopy  class 
representing  p  as  output. 


C(R(p)) 


ALG  2 


F(C(R(p))) 


F(R(p)) 


K(F(R(p))) 


K(F(R(p))) 


Figure  3.5   Flow  Chart  of  Two  Algorithms 

The  idea  of  the  test  is  to  determine  for  a  given  path 
its  raw  string  character  representation.  For  computational 
purposes,  the  paths  considered  are  polygonal  paths.  The  raw 
string  is  then  used  as  input  to  both  algorithms. 
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The  raw  string  R(p)  is  first  input  to  Algorithm  2. 
The  output  F(R(p))  from  Algorithm  2  is  then  input  to  the 
cancellation  routine  and  the  fully  cancelled  string  K(F(R(p))) 
is  output.  In  parallel,  R(p)  is  input  to  Algorithm  1  which 
produces  the  canonical  string  C(R(p)).  Algorithm  2  is  then 
applied  to  C(R(p))  to  output  F(C(R(p))).  The  results  of  these 
tests  are  then  compared.  In  Figure  3.5  a  flow  chart 
illustrates  this  process.  The  program  continues  to  compare 
the  two  outputs  until  a  predetermined  number  of  paths  are 
checked. 

It  can  be  shown  that  C(R(p))  is  a  unique  class  name  if 
and  only  if  F(C(R(p)))  equals  K(F(R(p)))  for  all  p  in 
T  [Ref .  2].  Before  the  proof  of  this  conjecture  was  obtained, 
the  above  procedures  were  programmed  and  successfully  tested 
for  one  million  different  paths. 
2 .   The  Test 

The  general  algorithm  used  to  test  the  two  procedures 
represented  in  Figure  3.5  is  easily  followed  and  is  shown  in 
Figure  3.6.  However,  the  attached  program,  which  was  used  to 
test  the  model,  is  substantially  more  involved.  This  more 
involved  program  was  written  in  an  attempt  to  generate  the 
most  efficient  programming  code  possible. 
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n  =  number  of  reference  frames  to  be  considered 
m  =  number  of  paths  to  be  considered 
nobs  =  number  of  obstacles  on  each  board 

begin 

for  board  =  1 ,  .  .  .  ,  n 

create  reference  frame  or  board 
for  path  =  1 ,  .  .  .  ,  n 

1.   Create  a  polygonal   path  p  with  nseg 

segments 
2  .   Form  R(p) 

3.  Form  F(R(p)) 

4.  Form  K(F(R(p))) 

5.  Form  C(R(p)) 

6.  Form  F(C(R(p))) 

if  K(F(R(p)))  does  not  equal  F(C(R(p)))  then 
print  (raw  string,  board,  and  path  info 
for  the  counter  example) 
end  if 
end  for 
end  for 

print  (final  seed) 
end 

Figure  3.6   Test  Algorithm 

Two  time  saving  techniques  are  employed  in  order  to 
test  the  one  million  cases  in  25  minutes.  The  first  method 
involves  the  subroutines  which  perform  cancellation.  It  is 
simple  to  code  a  program  that  scans  the  character 
representation  repeatedly  to  find  all  possible  cancellation  in 
a  string  such  as  B^B2ctiai^2/  wnicn  reduces  to  B^ .  However,  in 
the  enclosed  code  (Appendix  A),  pointers  are  inserted  in  the 
strings  to  mark  the  position  where  the  first  cancellation 
occurs.  Then  the  newly  adjacent  characters  are  checked  for 
possible  cancellation.  This  method  reduces  the  complexity  of 
the  test  and  is  employed  in  both  Algorithm  1  and  Algorithm  2. 
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The  second  time-saving  device  was  developed  around  the 
sorting  of  character  strings.  Instead  of  using  a  bubble  sort 
method  (order  n  complexity),  a  merge  sort  algorithm  (order 
n(log(n))  complexity)  is  used.  Table  1  shows  the  savings  is 
realized  immediately  as  opposed  to  there  existing  some  cut  off 
where  n(ln(n))  becomes  less  expensive. 


TABLE  1 
COMPLEXITY  COMPARISON 

n 

n2 

n(ln(n)) 

2 

5 

10 

100 

1,000 

1,000,000 

4 

25 

100 

10,000 

1,000,000 

10E16 

1 

8 

23 

460 

6,907 

13,815,510 

As  can  be  seen  in  the  attached  code,  this  choice  to 
save  computer  time  requires  a  significant  increase  in  program 
complexity  and  therefore  in  programmer  time.  The  trade-off  of 
computer  time  versus  programmer  time  needs  to  be  considered  in 
any  similar  analysis. 

The  code  used  to  implement  the  test  is  presented  in 
Appendix  A. 
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IV.   DETERMINING  A  BOUND  FOR  A  HOMOTOPY  CLASS 

A.   INTRODUCTION 

For  each  homotopy  class  on  a  list,  a  lower  bound  for  the 
length  of  its  shortest  representative  path  is  found.  The  list 
is  obtained  through  a  graph  traversal  process  [Ref.  3].  The 
list  of  classes  is  then  arranged  in  increasing  order  based 
upon  these  bounds. 

To  obtain  the  bounds,  we  first  fix  a  class  and  consider 
its  shortest  path,  p.  For  example,  consider  the  path  in 
Figure  4.1  to  be  a  string  attached  at  z.  The  shortest  path  is 
the  one  obtained  by  pulling  the  string  tight  so  that  it  lies 
against  the  obstacles  and  passes  through  "a". 


Figure  4.1   Class  Representatives 

The  topological  deformation  of  "contraction"  may  then  be 
applied  to  each  obstacle  to  "shrink"  the  obstacle  to  an  point 
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bjj  which  is  chosen  arbitrarily  within  the  obstacle  Bl  [Ref  .  4]. 
As  seen  in  Figure  4.2,  when  all  obstacles  can  be 
simultaneously  contracted  to  b^  .  The  shortest  path  p  found 
in  this  way  has  a  limiting  position  which  is  polygonal. 


Figure  4.2   Class  Representative  to  a  Bound  for 
the  Shortest  Representative 


The  length  of  this  polygonal  path  provides  a  lower  bound  that 
is  associated  with  that  class.  The  fact  that  this  is  a  lower 
bound  follows  from  the  following  argument.  If  we  let  D  be  the 
set  of  points  through  which  paths  may  travel  with  full-sized 
obstacles  and  E  be  the  set  of  available  points  in  the  space 
after  the  obstacles  are  "shrunk",  then  D  is  a  subset  of  E. 
Let  II(X,  E)  be  the  set  of  all  paths  in  equivalence  class  X 
which  remains  entirely  in  E.  Since  D  is  a  subset  of  E,  II(X,D) 
is  a  subset  of  II(X,E)  for  all  X.  We  now  define  a  function  f 
that  maps  H(X,E)  to  a  real  number  in  R,  with  the  rule  that 
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f(p)  =  length  of  p.   Since  D  is  a  subset  of  E,  this  equation 
foil ows : 


min  f(p)       min  f(p) 
peE  (X,  E)    ^  pell  (x,  D) 


(Figure  4.3) 


This  lower  bound,  min  pemj  ^  f  (p)  ,  can  be  calculated  from  the 
class  name  and  obstacle-representing  points  without  explicitly 
defining  any  contraction  deformation. 


n(X,E) 


Figure  4.3   Mapping  Related  to  Class  Lower  Bound 

Knowing  that  this  polygonal  path,  p  ,  exists,  we  now  give 
a  method  to  find  the  minimum  length  path  and  compute  its 
length.  Effectively,  we  must  connect  the  points  b^  in  the 
correct  order  and  compute  the  length  of  each  line  segment. 

The  idea  employed  to  find  p  is  :  We  are  given  a  starting 
point  "a"  and  a  character  string  which  represents  the 
mandatory  order  of  ray  crossings.  Each  character  of  the 
string  implies  that  we  must  cross  a  certain  ray.  We  then 
compute  the  maximum  number  of  rays  (i.e.  read  the  maximum 
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number  of  characters  in  the  string)  that  can  be  crossed  with 
a  straight  line.  If  z  cannot  be  reached  with  a  straight  line 
we  determine  where  it  is  best  to  make  a  turn.  The  turning 
point  then  becomes  the  new  starting  point  and  the  process 
repeats  until  z  is  reached.  To  find  the  straight  line 
segments  of  p  ,  the  notion  of  a  cone  of  directions  in  which 
the  segments  of  a  polygonal  p  must  lie  is  introduced. 

B.   FINDING  THE  BOUND 

1.   Establishing  a  Cone  of  Directions 

The  cone  of  directions  as  depicted  in  Figure  4.4,  is 
an  open  region--not  including  the  boundaries--f ormed  by  two 
rays  based  at  a  common  point  and  forming  an  angle  of  less  than 
%  radians.  In  this  region  we  can  draw  any  ray  based  at  the 
point  common  to  the  boundary  rays.  This  is  the  cone  of 
directions . 


Figure  4.4   Cone  of  Directions 


2.   Class  Name  Defines  the  Cone 

To  create  a  cone  we  need  a  starting  point  "a"  and  a 
reference  ray  (Figure  4.5).  We  define  an  entire  class  name  as 
Qi.  for  a  path  to  cross  the  ray  c^.   First  we  draw  a  ray  from 
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the  starting  point  to  the  point  representing  the  obstacle 
point  b^.  In  Figure  4.5,  to  cross  a^  we  must  be  to  the  right 
of  the  ray  from  a  to  b^.  The  ray  from  "a"  to  b^  is  one 
boundary  for  the  cone.  Next  we  draw  a  ray  parallel  to  the  ray 
that  we  must  cross,  o^ .  Again  in  Figure  4.5,  in  order  to 
cross  c^  with  a  straight  line  we  must  leave  "a"  at  an  angle 
that  will  allow  us  to  cross  the  Ql  ray.  Therefore,  the  ray 
parallel  to  Ql  forms  the  other  boundary.  The  two  boundaries 
form  the  cone  of  directions  as  any  ray  strictly  lying  between 
these  two  limits  will  cross  atL. 


starting 
point 


«k 


Figure  4.5   A  Cone  of  Directions  Toward  a^ 

3.   Intersecting  Multiple  Cones 

Given  a  cone  of  directions  for  a  single  ray,  we  can 
then  consider  longer  class  name  strings.  For  this  we  generate 
an  intersection  of  multiple  cones  of  direction.  For  example, 
in  Figure  4.6,  "a"  is  the  starting  point  and  the  class  is 
a2&\B2    leading  to  the  ending  point  z. 
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a2  _j 

Figure  4.6   The  Class  a^fli^ 


The  cone  of  directions  in  which  the  first  segment  of 
p  must  lie  is  determined  by  considering  the  first  element  of 
the  class  name.  Therefore,  the  first  cone  must  be  in  a 
direction  which  crosses  the  aj    reference  ray. 

We  draw  a  line  from  the  starting  point  a  to  bj  and  a 
second  line  from  "a"  parallel  to  c^  and  in  the  same  direction 
as  aj  (Figure  4.7).  This  forms  the  cone  of  directions  in 
which  the  first  polygonal  line  segment  of  the  shortest 
representative  must  lie. 


Figure  4.7   Cone  of  Directions  Toward  a2 
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Then  we  consider  the  second  character  B,     and  create 


1 
the  cone  of  directions  that  will  insure  crossing  of  B 


intersect   the   two   cones   of   direction, 
intersections  is  pictured  in  Figure  4.8. 


We 


The   cone   of 


new  cone 


intersection 


Figure  4.8   Two  Intersecting  Cones  of  Direction 

Continuing  in  this  fashion,  in  Figure  4.9,  we  consider  the 
third  character  in  the  class,  Q^  and  create  the  appropriate 
cone . 
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new  cone 


empty 
intersection 


Figure  4.9   An  Empty  Intersection  of  Directions 

We  see  that  the  intersection  of  these  cones  is  empty. 
To  enter  this  new  cone  we  are  forced  out  of  the  cone  in  which 
we  know  the  first  segment  must  lie.  Therefore,  we  choose  a 
path  which  crosses  all  rays  that  were  named  (a^  an<3  fii  )  and 
turns  on  some  bj..  We  determine  the  correct  b^  below. 
4.   An  Empty  Intersection 

As  can  be  seen  in  Figure  4.9,  the  new  cone  of 
directions  lies  to  the  right  of  the  previous  intersected  cone 
of  directions.  This  implies  that  the  path  we  require  will 
bend  to  the  right.  The  peg--point  b^  on  which  the  path  bends- 
-is  found  on  the  right  hand  boundary  of  the  old  cone.  In 
fact,  as  shown  in  Figure  4.9,  the  obstacle  centered  at  b^  is 
the  peg  that  the  path  bends  on  in  the  class  ctj^i^' 


31 


We  now  know  that  one  segment  of  the  bounding  path  is 
the  line  from  a  to  bj  .  In  the  next  step  we  use  bj  as  the  new 
starting  point.  Since  o^  and  Gj  have  been  crossed,  we  use  B2 
to  create  the  cone  of  directions  in  Figure  4.10. 


z  does  not 
lie  in  the  cone 


Figure  4.10   Cone  of  Directions 

Now,  since  there  are  no  more  characters  to  read,  we 
check  to  see  if  z  lies  within  the  most  recently  established 
cone  of  directions.  In  Figure  4.10  we  see  that  z  does  not  lie 
in  the  cone  so  the  bounding  path  must  bend  again. 

Since  z  lies  to  the  right  of  the  cone,  the  bounding 
path  bends  on  bj  which  is  on  the  right  hand  boundary  of  the 
cone  of  directions.  Hence,  the  final  two  segments  of  the 
bounding  path  are  the  line  segments  between  bj  and  b2  and  the 
line  segment  between  bj  and  z. 

So  given  the  class  ajBiB,  we  obtain  the  bounding  path 
pictured  in  Figure  4.11  on  the  contracted  obstacles. 
Computing  the  length  of  this  path  is  now  straight  forward. 
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The  next  section  generalizes  these  steps  and  shows  how 
to  find  the  path  to  z. 


Figure  4.11   The  Bounding  Path  for  the  Class  a2fllfl2 


THE  ALGORITHM 


To  find  a  lower  bound  for  a  path  in  a  given  class,  the 
following  algorithm  is  applied: 


t(j)     := 

peg(i) 
cone( t ( 

U    :  = 

Vold:  = 
V    :  = 

L,R    :  = 


Lrj,Rg 


an  array  of  points  with  t(l)  =  a,  and  t(m)  =  z; 
t(i+l)  through  t(m-l)  are  pegs  visited  along  the 
bounding  polygonal  path. 

:=    the  obstacle-representing  point  associated 
with  x^ 

j)/X:):=  the  cone  of  directions  from  a  starting 
point  t(j)  in  which  we  can  cross  the  ray 
x^  with  a  straight  line  from  t 
represents   the  cone  associated  with  the  most 
recently  read  character 

records  the  most  recent  cone  of  intersections 
represents  the  cone  of  directions  created  by  one 
character;  V  is  assigned  to  Vqj(j  and  is  replaced  by 
the  intersection  of  U  and  V  jj 

represent  pointers  into  the  character  string  X 
and  is  associated  with  the  obstacle  on  the  left 
of  right  boundary  of  the  cone  of  directions.  If 
no  boundary  lies  on  a  ray,  then  L  or  R  is 
assigned  zero 

associated  with  the  obstacle  on  the  left  of  right 
boundary  of  the  cone  U  that  is  intersected  with 


V 


old 


If  one  ray  of  U  lies  inside  V0u,    then  either 


L  or  R  is  replaced  by  Ln  or  R 


t- 


both  rays  of 
U  lie  inside  VoId  then  L'and  R  are  replaced  by  LD 
and   Rg    respectively. 
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BEGIN 

t(l)<-a;  i<-  1 
procedure  START 

read  character  \[ 

initialize  current  cone  V  <--  cone(t,  xj) 

IF  peg(  i )  on  left  boundary  of  V  THEN 

L<--i;    R<--0 
ELSE 

R<--i;    L<--0 
END  IF 
procedure  CONTINUE 
i  <--  i+1;  read  character  xj 

U  <--  cone(t,  xj) 

IF  peg(xj)  on  left  boundary  of  V  THEN 

Ltj<--i;    Rtj<--0 
ELSE 

Ru<--i;    Lu<--0 
END  IF 
Void  <"  V;  L0id  <--  L;  Roid  <--  R 

V  <-  U  n  Vold 

IF  left    boundary  of  U  falls  within  V0id  THEN;  L  <-- Ltj;  END  IF 

IF  right  boundary  of  U  falls  within  V0id  THEN;  R  <--  Rtj;  END  IF 

procedure  BRANCH 

IF  V  is  empty  THEN 

IF  U  is  to  left  of  V0ia  THEN 

j<--j  +  l;  t  (j)<-peg(L0id) 
ELSE  (  U  is  to  right  of  VGld  ) 

j<-j  +  l;  t  (j)<-  peg(R0id) 
END  IF 

go  to  procedure  START 
ELSE  (V  is  not  empty.) 

IF  i  <  m  (not  all  characters  have  been  read)  THEN  Go  to  procedure  CONTINUE. 
ELSE  (There  are  no  more  characters  to  read.) 

IF  z  lies  in  current  cone  V      THEN  j  <--  j  +  1;  t(j)  <--  z;  STOP 
ELSE 

IF  U  is  to  left  of  V0id  THEN 

j  <--  j  +  1;  t(j)  <--  peg(L0id);  j  <-  j  +  1;  t(j)  <--  z;  STOP 
ELSE  ( U  is  to  right  of  VGld  ) 

j  <--  j  + 1;  t(j)  <--  peg(Roid);  j  <--  j  + 1;  t(j)  <--  z;  stop 

END  IF 
END  IF 
END  IF 
END  IF 
END 
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V.   FINDING  THE  SHORTEST  PATH 

A.   INTRODUCTION 

In  order  to  find  the  true  shortest  path  from  a  to  z  we 
consider  the  classes  of  paths  on  the  ordered  list  of  candidate 
paths.  We  proceed  by  removing  the  first  class  from  the  list. 
The  true  shortest  path  in  the  named  class  is  found.  If  the 
length  of  this  path  is  smaller  than  the  bound  associated  with 
the  class  on  the  top  of  the  remaining  ordered  list,  the  search 
is  stopped.  Otherwise,  the  first  class  is  removed  from  the 
remaining  list  and  the  above  procedure  is  repeated.  This 
continues  until  the  true  shortest  path  is  found. 

To  find  the  shortest  path  within  a  given  class,  we  reverse 
the  contraction  applied  to  the  obstacles  performed  earlier 
when  we  found  the  lower  bound  for  the  classes  (see  Figure 
5.1).  For  each  polygonal  path,  whose  length  provided  a  lower 
bound,  we  transform  the  segment  into  the  true  shortest  path  in 
that  class.  Accordingly,  the  process  of  calculating  the 
shortest  path  in  a  class  begins  with  the  polygonal  path  used 
to  lower  bound  it  and  a  description  of  the  obstacle 
geometries . 


35 


Figure  5.1   Reversing  the  Obstacle  Contraction 

To  this  point,  no  assumptions  have  been  made  about  the 
shapes  of  the  obstacles.  In  this  chapter,  with  little  loss  of 
generality  and  for  computational  purposes  which  will  be 
discussed  later,  all  obstacles  are  assumed  to  be  circles 
centered  on  bj,  .  We  could,  of  course,  assume  other  shapes, 
like  polygons  or  ellipses. 

B.   FINDING  LINE  SEGMENTS 

As  can  be  seen  in  Figure  5.2  the  line  segment  from  a  to  B^ 
is  tangent  to  Bj.  Also,  the  line  segment  from  Bj  to  the 
obstacle  B2  is  tangent  to  both  obstacles.  In  general  there 
are  two  tangent  lines  from  a  point  to  a  circle  and  there  are 
four  distinct  lines  tangent  between  two  circles  (see  Figure 
5.2).  In  this  section  we  see  how  to  find  these  tangent  lines 
and  we  show  how  to  pick  the  correct  tangent  line  to  represent 
the  true  shortest  path  in  the  class. 
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Figure  5.2   Possible  Tangent  Lines 

Assuming  that  the  shortest  path  from  a  to  z  in  a  given 
class  of  paths  is  not  a  straight  line,  we  use  another 
application  of  the  cone  of  directions  idea  to  determine  on 
which  side  of  an  obstacle  to  turn.  That  is  equivalent  to 
picking  the  correct  tangent  line.  We  use  the  last  cone  of 
intersections  formed  just  before  the  intersection  of  cones 
becomes  empty  to  make  this  determination. 

In  Figure  5.3  we  see  that  the  path  is  required  to  bend 
around  obstacle  Bi  .  If  a  tangent  line  is  drawn  from  a  to  the 
left  side  of  obstacle  1  and  we  then  travel  on  an  arc  around 
the  obstacle  until  we  can  proceed  directly  to  z,  we  would  have 
the  shortest  path  in  the  class  and  could  easily  find  its 
length.   If  not  then  we  iterate  the  process. 
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Figure  5.3   Path  Bend  in  Class  B 


The  following  method  is  used  to  find  the  tangent  lines  to 
a  circular  obstacle  when  given  the  starting  point,  the 
coordinates  of  the  center  of  the  obstacle,  and  the  radius  of 
the  obstacle. 

In  Figure  5.4,  the  starting  point,  a,  is  the  point  (u,v). 
The  center  coordinates  and  the  radius  of  the  peg  are  the  point 
(h,k)  and  r,  respectively. 

y 


(h,k) 


>■   x 


a  (u,v) 

Figure  5.4   Point  to  Obstacle  Tangent  Lines 
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With  this  information  we  determine  the  equations  of  the 
circle  representing  the  obstacle  and  an  arbitrary  line.  We 
solve  these  equations  simultaneously  for  the  points  of 
intersection  (Figure  5.5).  The  specific  value  that  we  want  is 
the  point  of  tangency. 


Figure  5.5   The  Intersections  of  a  Circle  and  a  Line 


So,   given  h,  k,  and  r  and  assuming,  without  loss  of 
generality,  that  b=u=v=0,  we  find  the  following  two  equations: 

(x-h)'  +  (y-k)2  =  r2 

y  =  mx  (where  m  is  the  slope  of  the  line) 
By  expanding  the  first  equation  and  substituting  y=mx  we 
obtain : 

(l+m2)x2  +  (2mk-2h)x  +  (h2+k2-r2)  =  0  . 
By  the  quadratic  formula  we  obtain: 


x   = 


■2mk+2h±  [  (2mk-2h)  2-4  (1+n?)  (I?+k?-r2)  ]  * 
2(l+z#) 


(1) 


In  Equation  (1),  if  the  discriminant  is  zero,  we  have  found 
the  value  where  the  points  of  intersection  coincide.  So 
setting  the  discriminant  equal  to  zero  yields: 
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(-2A-2ion)2-4(l+i#)  (tf+tf'-r2)  =0 
(h+km)*-(nfi+l)  (tf+ic^-r2)  =0 
(k*n?+2hkm+&)  -(i^+l)  (&+k*-z2)  =0 
( ic'-tf-J^+r2)  iri*+  (2i2ic)  /n+  (r2-^)  =0 

from  which  m  can  be  determined. 

A  second  application  of  the  quadratic  formula  produces 

-2hk±  [4A2Jcg-4  (r*-W)  (r2-*2)  ]  * 
2(r2-A2) 

In   general ,    with   an   arbitrary    u   and   v 


m 


-  (h-u)  (k-v)  ±r  [  (h-u)  2+  (k-b)  2-r2]  *      ..    , 
i  "  r*-(h-u)*  '  ' 


The  only  time  this  technique  does  not  appear  to  work  is 
when  the  tangent  line  found  is  a  vertical  line  (m  is  undefined 
because  there  is  a  zero  in  the  denominator).  This  problem  can 
be  easily  handled  because  if  m  is  undefined  we  know  the 
equation  for  a  vertical  line  to  be  x  equal  to  some  constant. 

The  value  for  m  is  then  substituted  into  Equation  (1)  to 
find  the  xl ,  i  =  1,2,  coordinate  and  the  value  of  Xj  is 
substituted  into  y^  =  rn^X'  to  find  y^.  By  knowing  on  which  side 
of  the  obstacle  the  path  must  turn,  we  can  pick  the  tangent 
line  with  the  appropriate  slope.  A  method  of  picking  the 
correct  line  will  be  discussed  in  the  next  section. 

40 


We  follow  the  arc  along  the  circle  to  the  point  of 
tangency  of  the  line  between  the  current  obstacle  and  z  or  to 
the  point  of  tangency  to  the  next  obstacle.  Tangent  lines 
between  two  obstacles  can  be  determined  by  a  similar  method  as 
the  one  just  shown.  A  method  of  finding  the  four  tangent 
lines  between  two  circular  obstacles  follows. 


C.   DETERMINING  THE  CORRECT  TANGENT  LINES 
1.   Point  to  Obstacle 

We  have  shown  how  to  find  the  two  tangent  lines  from 
a  point  to  a  circle.  We  present  a  method  to  determine  which 
tangent  line  to  pick. 

The  obstacle  on  which  we  turn  lies  on  a  boundary  of 
the  cone  of  directions.  Thus  the  cone  itself  lies  either  to 
the  left  or  to  the  right  of  the  ray  through  the  obstacle 
center.  If  the  cone  lies  on  the  left,  use  the  tangent  of  the 
left  is  taken.  If  the  cone  lies  to  the  right,  then  the 
tangent  on  the  right  is  taken.  The  case  with  the  cone  on  the 
left  is  illustrated  in  Figure  5.6. 


Figure  5.6   Use  the  Left  Tangent  for  Path  Bending  to  the  Right 
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2.   Obstacle  to  Obstacle 

When  considering  two  circular  obstacles,  there  are 
four  distinct  lines  that  can  be  drawn  tangent  to  both 
obstacles.  In  this  section  we  show  how  to  pick  the  tangent 
line  that  corresponds  to  the  respective  shortest  path  in  the 
given  cl ass . 

To  find  the  four  distinct  tangent  lines  between  two 
obstacles  we  first  create  a  line  which  is  tangent  to  one  of 
the  circles.  We  then  move  this  tangent  line  along  the  circle 
boundary  until  we  discover  the  four  points  at  which  the  line 
is  tangent  to  both  obstacles  (Figure  5.7). 


(x',y")  =  (r  cos(<|>),  r  sin(<») 


->-  x 


y  =  mx  +  b 


Figure  5.7   Rotating  $   to  Find  Tangent  Lines  Between 
Two  Circles 
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Without  loss  of  generality,  in  Figure  5.7  we  fix  one 
obstacle  center  at  the  origin.  We  pick  an  arbitrary  point 
(x',y'),  on  the  obstacle  Bj  represented  as  (r*cos(<j>), 
r*sin(<|>)).  We  construct  a  line  through  (x',y')  tangent  to  Bi  . 
The  slope  m  of  this  line  is: 

y-r*sin(<|>) 
x-r*cos  (<J>) 

This  slope  is  also  equal  to  -cot(<|>).   Hence, 

m=   y-r*sin<|>   =  _cot  (<())  (2) 

x-r*cos<|> 

The  equation  of  the  circle  with  center  (h,k)  is 
(x-h)2+(y+k)2   =  R2 

Rearranging  Equation  (2)  we  obtain  the  following: 

y  -  r*sin(<|>)  =  (-cot  (4>)  )  (x  -  r*  (cos  (<t>) )        ,  v 

(3) 

y=  (-cot  (<|))  )x  +  r*cos  (<|>)  *cot  (<|>)  +  r*sin(<|>) 

Equation  (3)  is  in  slope-intercept  form  and  the  y  intercept  b 
is  : 
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=  r 


b  =  r*cos  (<|>)  *cot  (<|>)    +  r*sin(4>) 

[S2#±  +Sin*)) 

^  sm<f>  / 

/  cos2^  +  sinty  \ 
{  sin<|>  J 


=  r 
b  =  r*csc(<|>) 

Thus  both  m  and  b  are  functions  of  <(>.  The  general  line  that 
we  use  to  intersect  both  circles  is 

y  =  m($)x  +  i>(4>)  . 

This  equation  for  the  line  and  for  the  obstacle  centered  at 
(h,k)  are  solved  simultaneously.  After  steps  similar  to  those 
used  in  the  point  to  obstacle  example  we  obtain  the  following 
formula: 

(/n*+l)  x*   +  [-2h+2m(b-k)]x  +  [#*+  (jb-ic)  *-&]    =0  . 

Again  by  applying  the  quadratic  formula  and  setting  the 
discriminant  equal  to  zero  we  find 

4[-h-m(b-k)]*  -  4(2n*+l)  [#+ (jb-ic)  *-;?]  =0  . 

Substituting  the  function  values  for  m(<|>)  and  b(<|>)  yields 
[-h-(cotQ)  (r*csci>-k)  ]*-(cot*$+l)  [fr+(r*csci>-k)*-&]  =0  (4) 
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This  equation  has  four  zeros  for  values  of  (j)  in  [0,2x) 
as  can  be  seen  in  Figure  5.8. 
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Figure  5.8   Four  Zeros  of  Equation  (4) 

Since  <|>  is  the  only  parameter  in  Equation  (4)  we  can 
use  the  four  values  found  to  label  the  four  tangent  lines. 
The  line  associated  with  the  lowest  value  of  $  is  labeled  A, 
the  second  line  associated  with  the  second  lowest  4>  value  is 
labeled  B  and  the  last  two  lines  are  labeled  with  C  and  D, 
respecti vel y . 
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Y  c 

Figure  5.9   Possible  Tangent  Lines 


>-  X 


Assuming  the  polygonal  path  bends  on  each  peg  on  the 
path,  there  are  essentially  two  possible  arrangements  when 
considering  a  polygonal  path  that  negotiates  two  obstacles. 
These  possibilities  are  shown  in  Figure  5.10.  The  paths  can 
be  flipped  about  the  horizontal  axis  to  produce  two  equivalent 
paths . 


Figure  5.10   Two  Possible  Polygonal  Paths 
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In  order  to  discover  which  situation  exists,  we  extend 
the  middle  line  segment  and  use  it  as  a  reference  in  Figure 
5.10.  Since  the  path  bends  at  each  peg,  the  two  end  segments 
must  lie  on  one  side  or  the  other  of  the  reference  line.  If 
they  both  lie  on  the  same  side,  then  we  have  situation  a 
above.  If  the  two  segments  lie  on  opposite  sides,  then  we 
have  situation  b. 

Given  the  two  possible  arrangements  in  Figure  5.10,  we 
know  the  corresponding  shortest  paths  in  the  class  look  like 
those  in  Figure  5.11. 


Figure  5.11   True  Shortest  Paths 
We  have  shown  how  to  pick  the  correct  tangent  line 
from  a  point  to  a  circle.   This  is  applied  to  find  the  first 
and  last  segments  of  this  polygonal  path. 

To  find  the  correct  tangent  between  two  obstacles, 
following  the  situation  on  the  left  in  Figure  5.11  we  must 
pick  tangent  line  B,  described  earlier.  If  we  have  the  other 
situation  in  Figure  5.11,  we  must  pick  tangent  line  A. 

If  the  polygonal  path  does  not  bend  on  each  peg,  then 
there  are  two  possible  configurations  as  seen  in  Figure  5.12. 
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Figure  5.12   Polygonal  Paths  With  Collinear  Segments 


To  this  point,  we  only  used  the  cone  of  directions 
idea  to  find  the  correct  tangent  from  a  point  to  an  obstacle. 
We  use  the  cone  of  directions  information  to  pick  the  correct 
tangent  lines  too. 

As  seen  in  Figure  5.13  we  consider  the  class  B^aj  . 


—  bound 


true 


~1 

Figure  5.13   Class  BjCtj  and  its  Bounds 


The  tangent  line  from  a  to  Bi  is  found  in  the  same 
manner  as  discussed  earlier.  However,  since  the  path  of  the 
bound  does  not  bend,  we  cannot  choose  the  correct  tangent  line 
between  Bi  and  B*  in  the  same  way  as  before. 

Instead,  we  construct  a  reference  line  that  connects 
the  centers  of  the  two  obstacles.  Then  we  consider  on  which 
side  of  this  line  the  cones  of  directions  lie  (Figure  5.14). 
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If  they  lie  on  different  sides  of  the  reference  line  then  we 
pick  tangent  line  A  or  D.  Otherwise,  we  pick  tangent  lines  B 
or  C . 


empty 
intersection 


Figure  5.14   Cone  of  Direction  to  Help  Pick  Tangent  Lines 

D.   TWO  POTENTIAL  PROBLEMS  TO  CHECK 

Given  a  path  in  a  class  like  the  one  in  Figure  5.15,  we 
see  the  associated  bounding  path  when  we  shrink  the  obstacles 
to  points.  However,  there  are  two  separate  problems  that  may 
arise  when  the  obstacles  are  restored  to  their  original  shape: 
1)  The  prospective  path  may  pass  through  one  or  more 
obstacles.  2)  When  the  obstacles  are  restored  to  their  full 
size  the  path  created  can  switch  sides  of  one  or  more 
obstacles.   These  possibilities  can  be  seen  in  Figure  5.16. 
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Figure  5.15   A  Class  and  Its  Bounding  Path 
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the  bounding 
path 


Figure  5.16   A  Bounding  Path  and  Its  Restored  Path 

The  path  determined  by  the  cone  of  directions  technique 
passes  through  obstacles  B,  and  Bj,  illustrating  the  first 
problem.  By  switching  sides  of  B*,  the  path  is  no  longer  in 
the  same  equivalence  class,  illustrating  the  second  problem. 
Actually,  the  path  should  look  like  the  one  in  Figure  5.17 
which  avoids  all  obstacles  and  remains  in  the  original  class. 
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Figure  5.17   Corrected  Shortest  Path 

The  correct  path  around  the  obstacles  can  be  determined  by 
using  a  distance  function.  Then,  methods  discussed  in  Chapter 
IV  are  employed  to  negotiate  the  new  impeding  obstacles.  In 
Figure  5.18,  segment  Ps  --which  is  one  segment  of  the  bounding 
path--joins  pegs  b^  and  bwj,  and  lies  along  the  line  through 
these  pegs.   We  let  a  circular  obstacle  Bj  have  center  (u-,v) 
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and  radius  r„-  .  Pegs  b^  and  b^  have  cartesian  coordinates 
(x;,y:)  and  (xi  +  i'yi  +  l^  respectively.  We  let  d(*,*)  denote 
Euclidean  distance. 

The  distance  from  p   to  B>--which  is  the  perpendicular 

-1       J 

distance  from  ps  to  the  nearest  point  in  B-p-is  given  by 
d(p  ,Bj)  if  the  center  of  B^  lies  in  the  band  denoted  in  Figure 
5.18  whose  borders  are  perpendicular  to  p  and  which  passes 
through  b^  and  b^,  ,  and  by  min{d(bj[ ,  E-)  ,  d(  b^  ,  B  •)  }  if  not.  We 
can  then  check  whether  the  center  of  B.  is  in  the  band  if  we 


define  vectors 


and 


xi  =  (C*i+r*i]'  [yiii-yi]) 


u^    =    ([u-x^,    [vj-yj]) 


which  can  be  seen  in  Figure  5.19  [Ref.  5] 


not 
violated 


Figure  5.18   Band  Created  By  One  Segment  of  a  Polygonal  Path 
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Then  the  center  of  B]  is  in  the  band  if  and  only  if  X:  •  u,_i  ^>0 
and  X-    •    p-<0,  where  (•)  denotes  the  inner  product.   We  write 

the  distance  function  corresponding  to  p   passing  through  Bj 

5  J 

as 


dipl.Bj)     = 


dip^bj)  ;     if    \±  •  \i±_ltj>0       and    X±  •  ni:/<0 
min  d(bk,Bj)  ,  ^(£^,.8.7)  ;  otherwise 


These    distances    are    easily    computed   with 


d{bk.x.B3)    =    [(xruj)^  +    {yt-Vj)*]"  -  rj 


and 


dipl.Bj)- 


Kyitl-yi)yi  -    (x^-x^Vj-  x^u^-h^  +  y^x^-x^ 


-  i. 


[Ref.     5]. 


Figure  5.19   Vectors  to  Help  Determine  d(ps  ,Bj) 
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With  this  new  information,  appropriate  tangent  lines  must 
be  found  in  order  to  avoid  newly-violated  obstacles  and  to 
keep  the  path  in  the  correct  class.  These  new  tangent  lines 
may  also  encroach  obstacles.  So  the  process  of  checking  all 
obstacles  is  continued  until  no  obstacles  are  violated. 

To  find  the  length  of  the  shortest  path,  we  sum  the 
straight  line  segments  and  all  of  the  arc  lengths 
corresponding  to  the  distances  traveled  around  each  circular 
obstacl e . 
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VI.   CONCLUSIONS  AND  RECOMMENDATIONS 

A.   CONCLUSIONS 

Determination  of  the  shortest  path  between  two  points  in 
a  plane  containing  obstacles  has  applications  in  many  fields, 
particularly  that  of  robotic  path  planning.  This  thesis 
addresses  several  key  issues  belonging  to  a  larger  plan  for 
the  solution  of  the  shortest  path  problem.  This  plan  may  be 
summarized  as  follows: 

1.  Partition  the  set  of  paths  into  homotopy  classes. 

2.  Associate  algebraic  "names"  with  each  of  these  classes. 

3.  Produce  a  finite  list  of  feasible  classes. 

4.  For  each  class  on  the  above  list,  calculate  an 
inexpensive  lower  bound  on  the  length  of  its  shortest 
path. 

5.  Arrange  the  list  in  increasing  order  of  these  lower 
bounds . 

6.  Considering  the  classes  in  listed  order,  calculate  the 
shortest  path  in  the  class.   Calculation  terminates  when 
some  class  is  encountered  whose  shortest  path  has  length 
shorter  than  the  lower  bound  associated  with  the 
succeeding  class. 

The   approach   outlined   above   offers   two   primary 

computational  advantages.   First,  note  that  the  collection  of 

possible  paths  from  a  given  point  of  origin  to  a  given 

destination  is  uncountably  infinite.  By  organizing  the  search 

around  homotopy  classes  and  by  taking  advantage  of  certain 

topological  relationships  in  the  region  [Ref.  3],  the  search 

of  an  uncountably  large  set  is  replaced  by  a  search  of  a 

finite  list . 
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Second,  the  search  of  this  finite  list  is  ordered 
according  to  lower  bounds  which  estimate  the  lengths  of  the 
associated  shortest  paths.  This  ordering  allows  the  search  to 
terminate  without  exhausting  the  list.  The  bounds  associated 
with  each  class  are  inexpensive  to  calculate—relative  to  the 
cost  of  finding  the  class'  shortest  path.  This  results  in 
further  savings  of  computational  effort. 

The  contributions  of  this  thesis  concentrate  in  areas  2, 
4,  and  6  of  the  above  list. 

As  regards  the  association  of  names  with  homotopy  classes 
(area  2)  we  present  a  computational  investigation  which  was 
conducted  ancillary  to  proof  activity  in  this  area.  In 
summary,  we  associate  a  name  with  each  homotopy  class  as 
follows:  by  imposing  a  certain  labeled  reference  frame  on  the 
region  in  question;  by  encoding  information  about  the 
relationship  of  a  given  path  with  the  frame;  and  by  employing 
algorithms  presented  here.  In  this  way,  we  reduce  the  encoded 
information  to  a  string  of  characters  which  is  unique  to  each 
class.  The  investigation  presented  here  uses  a  presentation 
of  the  well-known  fundamental  group  as  a  standard  of 
comparison  to  verify  that  the  character  strings  used  in  our 
approach  are  in  fact  class  names.  (This  proposition  has 
subsequently  been  proved  [Ref.  2]).  The  significance  of  the 
structure  of  the  names  which  we  employ  is  that  they 
incorporate   information   which   is   useful   both   for   the 
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calculation  of  class  bounds  and  for  calculation  of  class 
shortest  paths. 

For  the  calculation  of  class  bounds  (area  4)  we  represent 
each  obstacle  by  a  point  chosen  within  it.  This  point  is  used 
to  construct  the  reference  frame  which  is  imposed  on  the 
region.  Restricting  our  attention  to  such  summary  information 
allows  class  bounds  to  be  calculated  cheaply. 

The  final  step  (area  6)  begins  with  an  ordered  list  of 
classes.  In  this  paper  we  present  a  method  by  which  shortest 
paths  of  a  class  may  be  calculated  when  the  obstacles  are 
circular.  This  is  the  only  step  in  which  any  generality  is 
lost,  and  is  potentially  fruitful  subject  for  further  study. 

B.   RECOMMENDATIONS  FOR  FURTHER  STUDY 

The  field  of  robotics  is  one  that  is  gaining  much  interest 
in  mathematics  and  other  fields  of  study.  Further  analysis  of 
the  reference  frame  and  its  properties  may  aid  future 
researchers  in  answering  other  questions  concerning  the 
shortest  path  problem. 

One  aspect  of  the  problem  that  can  be  studied  further 
relates  to  areas  4  and  6  above.  Programming  the  two  steps  of 
finding  the  lower  bound  for  a  class  and  subsequently 
determining  the  shortest  path  in  a  class  would  be  quite 
useful . 

Another  computational  problem  relates  to  area  2.  It  is 
believed  that  further  study  and  refinement  of  the  code  could 
result  in  a  program  which  is  even  more  efficient. 

58 


Another  possible  subject  of  future  research  is  that  of 
considering  non-circular  obstacles  (area  6).  It  should  be 
noted  that  any  set  of  obstacles  can  be  covered  with  circles. 
With  this  approximation  the  procedures  presented  in  this 
thesis  may  be  applied. 

Finally,  this  thesis  introduces  only  part  of  the  solution 
to  the  shortest  path  problem.  The  thesis  by  CAPT  Kevin  D. 
Jenkins,  U.S.  Marine  Corps  [Ref.  3],  presents  the  remaining 
portions  of  the  solution  (area  3),  and  should  be  studied  in 
conjunction  with  this  paper. 
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APPENDIX  A.   FORTRAN  PROGRAM:  THE  COMPUTATIONAL  INVESTIGATION 


PROGRAM  CLASNAME  FORTRAN 

C  THIS  PROGRAM  RANDOMLY  GENERATES  OBSTACLES  AND  A  POLYGONAL  PATH 

C  THROUGH  THOSE  OBSTACLES  TO  DETERMINE  WHETHER  OR  NOT  THE  PROCEDURES 

C  USED  BY  COMPETING  ALGORITHMS  PRODUCE  THE  SAME  FUNDAMENTAL  GROUP 

C  REPRESENTATION  FOR  THAT  PATH. 

C  THIS  PORTION  OF  THE  PROGRAM  SERVES  AS  THE  MAIN  DRIVER  WHICH  RECEIVES 
C  THE  PARAMETERS  DEFINING  THE  REGION  AND  INITIALIZES  THE  LINK  LIST 
C  ARRAYS  WHICH  WILL  BE  USED  TO  REPRESENT  THE  POLYGONAL  PATH. 

C  INPUT:   INITIAL  SEED  FOR  THE  RANDOM  NUMBER  GENERATOR,  NUMBER  OF 
C         OBSTACLE  CONFIGURATIONS  AND  PATHS  TO  BE  TESTED,  NUMBER  OF 
C         OBSTACLES  IN  THE  PLANE  AND  NUMBER  OF  SEGMENTS  IN  EACH 
C         POLYGONAL  PATH 

C  OUTPUT:   FINAL  RANDOM  NUMBER  GENERATOR  SEED  AND  MESSAGES  INDICATING 
C  ANY  PATHS  WHICH  PRODUCE  DIFFERENT  FUNDAMENTAL  GROUP 

C  REPRESENTATIONS 

REAL*8   BX(1000),  BY(1000),  X(1000),  Y(1000),  DSEED 
INTEGER  NOBS,  NSEGS,  NUMPTS,  N,  M,  HEAD(1000) 
+        NEXT (1000),  PRED(1000),  BOARD,  PATH 
EXTERNAL  GGUBFS 

DSEED  =  123457.0 

M  =  1000 

N  =  1000 

NOBS   =  20 

NSEGS  =  5 

NUMPTS  =  NSEGS  +  1 

PR  I  NT* , ' I NPUT  SEED ' ,  DSEED 

DO  1  BOARD  =  1,  N 

CALL  BOARDS  (DSEED,  BX,  BY,  NOBS) 
DO  2  PATH  =  1,  M 

CALL  IN  IT  (HEAD,  NEXT,  PRED,  PSEED,  DSEED  ) 
CALL  PATHS  (DSEED,  X,  Y,  NSEGS) 

CALL  TEST ( HEAD , NEXT , PRED , BX , BY , NOBS , X , Y , NUMPTS , BSEED , PSEED ) 
2     OONT I NUE 
1  CONTINUE 

PRINT*,*  ' 

PR  I  NT*,' FINAL  SEED',  DSEED 

STOP 
END 
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SUBROUTINE  I  NIT  (  HEAD,  NEXT,  PRED  ,  PSEED,  DSEED) 

C  THIS  SUBROUTINE  INITIALIZES  THE  DOUBLE  LINK  LIST  ARRAYS  THAT  WILL 
C  REPRESENT  THE  PATH. 

C   INPUT:   DSEED 

C  OUTPUT:  HEAD,  NEXT,  AND  PRED  ARRAYS  SET  TO  ZERO,  AND  SEED  FOR  THE 
C  R.N.G.  PRIOR  TO  CONSTRUCTION  OF  THE  PATH 

INTEGER  HEAD(1000),  NEXT(1000),  PRED(1000) 
REAL*8  PSEED,  DSEED 

PSEED  =  DSEED 

DO  3  K  =  1,1000 
HEAD(K)  =  0 
NEXT(K)  =  0 
PRED(K)  =  0 
3  CONTINUE 

RETURN 
END 


SUBROUTINE  BOARDS  (DSEED,  BX,  BY,  NOBS) 

C  THIS  SUBROUTINE  USES  A  PSEUDO  RANDOM  NUMBER  GENERATOR  TO  CREATE  THE 
C  COORDINATES  OF  EACH  OBSTACLE  ON  THE  BOARD,  SCALING  ALL  COORDINATES 
C  TO  BE  IN  THE  INTERVAL  (-1,1). 

C   INPUT:   NUMBER  OF  OBSTACLES  IN  THE  REGION  AND  A  SEED  FOR  THE  R.N.G. 

C  OUTPUT:   OBSTACLE  COORDINATES 

REAL*8  BX(1000),  BY (1000),  DSEED 
EXTERNAL  GGUBFS 

DO  1  I  =1,  NOBS 

BX(I)  =  2.  *  GGUBFS(DSEED)  -  1. 
1     BY(I)  =  2.  *  GGUBFS(DSEED)  -  1. 

RETURN 
END 
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SUBROUTINE  PATHS  (DSEED,  X,  Y,  NSEGS) 

C  RANDOMLY  GENERATES  THE  X  AND  Y  COORDINATES  OF  THE  VERTICES  FOR  THE 

C  POLYGONAL  PATH,  SCALING  ALL  COORDINATES  TO  BE  IN  THE  INTERVAL  (-1,1). 

C  THIS  SUBROUTINE  ALSO  ENSURES  THAT  THE  PATH  IS  A  CLOSED  LOOP  BY 

C  ASSIGNING  THE  START/FINISH  POINTS  THE  SAME  COORDINATES. 

C   INPUT:   NUMBER  OF  PATH  SEGMENTS  AND  A  SEED  FOR  THE  R.N.G. 

C  OUTPUT:   COORDINATES  OF  VERTICES  ALONG  THE  POLYGONAL  PATH 

REAL*8  X(1000),  Y(1000),  DSEED 
EXTERNAL  GGUBFS 

DO  1  I  =1,  NSEGS 

X(l)  =  2.  *  GGUBFS(DSEED)  -  1. 
1     Y(l)  =  2.  *  GGUBFS(DSEED)  -  1. 

X(NSEGS+1)  =  X(1) 
Y(NSEGS+1)  =  Y(1) 

RETURN 
END 


SUBROUTINE  TEST (HEAD,  NEXT,  PRED,  BX,  BY,  NOBS,  X,  Y,  NUMPTS, 

+  BSEED,  PSEED) 

C  ************************************************************ 

C  GENERATES  THE  RAW  STRING  OF  CHARACTERS  REPRESENTING  A  PATH  AND 

C  DETERMINES  ITS  FUNDAMENTAL  GROUP  REPRESENTATION  USING  ALGORITHMS  1 

C  AND  2.   RESULTS  OF  THESE  COMPETING  ALGORITHMS  ARE  THEN  COMPARED  FOR 

C  DIFFERENCES. 

C   INPUT:   A  POLYGONAL  PATH  AND  COORDINATES  OF  POINTS  REPRESENTING 
C         OBSTACLES 

C  OUTPUT:   COMPARED  RESULTS  OF  ALGORITHMS 

INTEGER  HEAD(1000),  NEXT(1000),  PRED(1000),  FR(1000),  NR(1000), 
+        PR(1000),  FL(1000),  NL(1000),  PL(1000) 
REAL*8  BX(1000),  BY(1000),  X(1000),  Y(1000) 

DATA  FR/ 1000*0/ 
DATA  NR/ 1000*0/ 
DATA  PR/ 1000*0/ 
DATA  FL/ 1000*0/ 
DATA  NL/ 1000*0/ 
DATA  PL/ 1000*0/ 
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CALL  RAWSTR(BX,  BY,  NOBS,  X,  Y,  NUMPTS,  HEAD,  NEXT,  PRED,  NELEMS) 

CALL  ALG2  (NOBS,  NELEMS,  HEAD,  NEXT,  X,  Y,  BX,  BY,  FR,  NR,  PR) 

CALL  CANALG2  (FR,  NR,  PR) 

CALL  ALG1  (HEAD,  NEXT,  PRED,  NELEMS) 

CALL  ALG2  (NOBS,  NELEMS,  HEAD,  NEXT,  X,  Y,  BX,  BY,  FL,  NL,  PL) 

CALL  CTREX  (FL,  NL,  PL,  FR,  NR,  PR,  BSEED,  PSEED) 

RETURN 
END 


SUBROUTINE  RAWSTR(BX,  BY,  NOBS,  X,  Y,  NUMPTS,  HEAD,  NEXT,  PRED, 
+  NELEMS) 

C  PRODUCES  THE  RAW  STRING  OF  CHARACTERS  WHICH  REPRESENTS  THE  PATH, 

C  PAYING  CLOSE  ATTENTION  TO  THE  ORDERING  OF  THE  CHARACTERS  WHERE 

C  NECESSARY.   THE  STRING  IS  CONSTRUCTED  BY  IDENTIFYING  THOSE  OBSTACLE 

C  REFERENCE  LINES  WHICH  ARE  CROSSED  AS  EACH  PATH  SEGMENT  IS  TRAVERSED 

C  IN  ORDER. 

C   INPUT:   ORDERED  LIST  OF  VERTICES  REPRESENTING  THE  PATH  AND  A  POINT 
C         REPRESENTING  EACH  OBSTACLE  IN  THE  REGION 

C  OUTPUT:   A  RAW  STRING  OF  ALPHAS  AND  BETAS  CONTAINED  IN  AN  ARRAY 

C  NAMED  HEADQ  AND  ITS  PARALLEL  ARRAYS  NEXT()  AND  PRED()  WHICH 

C  PRODUCE  THE  DOUBLE  LINKED  LIST 

IMPLICIT  REAL*8  (A-H,0-Z) 

INTEGER  HEAD,  NEXT,  SEG,  FSTSTR,  HDINDX,  SEGEND,  FIRST,  START, 
+        LENGTH,  PRED 

LOGICAL  ALLALF 

DIMENSION  X(1000),  Y(1000),  BX(1000),  BY(1000),  A1(1000), 
+         B1(1000),  A2(1000),  B2(1000),  D2(1000),  DIST(1000), 
+         FIRST(IOOO),  HEAD(1000),  NEXT(1000),  PRED(1000) 

DATA  D I  ST/ 1000*0.0/ 
DATA  FIRST/ 1000*0/ 

CALL  SETUP(NUMPTS,NOBS,CX,CY,X,Y,BX,BY,A1 ,B1 ,A2,B2,D2) 

HDINDX  =  2 
NSEGS  =  NUMPTS  -  1 
LOLD  =  1 
HEAD(LOLD)  =  0 
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DO  7  SEG  =  1,  NSEGS 

SEGEND  =  SEG  +  1 

FSTSTR  =  HDINDX 

LSTSTR  =  HDINDX 

ALLALF  =  .TRUE. 

LENGTH  =  0 

DO  6  LINE  =  1,  NOBS 

CHECK1  =  (A1(LINE)*X(SEG)+B1(LINE)*Y(SEG))* 

+  (A1 (L I NE )*X( SEGEND )+B1 (L I NE)*Y(SEGEND) ) 

IF  (CHECK1.LT. 0)  THEN 
LENGTH  =  LENGTH  +  1 
CHECK2  =  (A2(SEG)*BX(LINE)+B2(SEG)* 
+  BY(LINE)+D2(SEG))*D2(SEG) 

IF  (CHECK2.LT. 0)  THEN 

HEAD(HDINDX)  =  -LINE 
LSTSTR  =  HDINDX 
HDINDX  =  HDINDX  +  1 
ELSE 

CALL  CASES1  (A1,  B1,  A2,  B2,  D2,  SEG,  LINE,  XINT, 
+  YINT) 

DISTC  =  XINT**2  +  YINT**2 

DISTB  =  (XINT-BX(LINE))**2  +  (YINT-BY(L INE))**2 
IF  (DISTB. LT. DISTC)  THEN 
HEAD(HDINDX)  =  LINE 
LSTSTR  =  HDINDX 
HDINDX  =  HDINDX  +  1 
ALLALF  =  .FALSE. 
ELSE 

HEAD(HDINDX)  =  -LINE 
LSTSTR  =  HDINDX 
HDINDX  =  HDINDX  +  1 
END  IF 
END  IF 
END  IF 

6  CONT I NUE 

IF  (LENGTH. NE.O)  THEN 
IF  (ALLALF)  THEN 

CALL  ALPHAS  (NEXT,  LOLD,  FSTSTR,  LSTSTR,  SEG,  START, 
+  HEAD,  NSEGS,  PRED) 

ELSE 

CALL  ORDER(SEG,  LINE,  HDINDX,  HEAD,  NEXT,  A1,  A2, 
+  B1,  B2,  D2,  FSTSTR,  LSTSTR,  X,  Y,  START, 

+  LENGTH,  FIRST,  LOLD,  NSEGS,  PRED) 

END  IF 
END  IF 

7  CONTINUE 

CALL  COUNTR( START,  NEXT,  HEAD,  NELEMS) 

RETURN 

END 
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SUBROUTINE  SETUP ( NUMPTS ,  NOBS,  CX,  CY,  X,  Y,  BX,  BY,  A1,  B1,  A2, 

+  B2,  D2) 

C  ********************************************************* 

C  FOR  EACH  OBSTACLE,  THIS  ROUTINE  DETERMINES  THE  COEFFICIENTS  OF  THE 

C  EQUATION  FOR  THE  REFERENCE  LINE  FROM  THE  OBSTACLE  TO  THE  ORIGIN.   IN 

C  ADDITION,  IT  CALCULATES  THE  COEFFICIENTS  OF  THE  LINE  REPRESENTING 

C  EACH  SEGMENT  OF  THE  POLYGONAL  PATH. 

C   INPUT:   NUMPTS,  NOBS,  COORDINATES  OF  VERTICES  ALONG  POLYGONAL  PATH 
C         AND  COORDINATES  OF  THOSE  POINTS  WHICH  REPRESENT  EACH 
C         OBSTACLE 

C  OUTPUT:   COEFFICIENTS  OF  LINEAR  EQUATIONS  REPRESENTING  PATH  SEGMENTS 
C  AND  REFERENCE  LINES  FOR  EACH  OBSTACLE 

IMPLICIT  REAL*8  (A-H,  O-Z) 

DIMENSION  X(1000),  Y(1000),  BX(1000),  BY(1000),  A1(1000), 
+         B1(1000),  A2(1000),  B2(1000),  D2(1000) 

DO  4  I  =  1,  NOBS 
A1(l)  =  BY(I) 
B1(l)  =  -BX(I) 

4  CONTINUE 

CX  =  0.0 
CY  =  0.0 

I  =  1 

DO  5  J  =  2, NUMPTS 

A2(l)  =  Y(J)  -  Y(l) 

B2(l)  =  X(l)  -  X(J) 

D2(l)  =  (Y(I)*X(J))-(X(I)*Y(J)) 

I  =  1  +  1 

5  CONTINUE 

RETURN 

END 
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SUBROUTINE  CASES1(A1,  B1,  A2,  B2,  D2,  SEG,  LINE,  XINT,  YINT) 

C  THIS  SUBROUTINE  DETERMINES  THE  COORDINATES  FOR  THE  POINT  OF 

C  INTERSECTION  OF  A  GIVEN  PATH  SEGMENT  AND  A  GIVEN  OBSTACLE  REFERENCE 

C  LINE.   NUMERICAL  STABILITY  OF  CALCULATIONS  REQUIRES  THE  MATHEMATICAL 

C  OPERATIONS  TO  BE  SEPARATED  INTO  CASES,  THE  BEST  CASE  BEING  USED  FOR 

C  EACH  PARTICULAR  SITUATION. 

C   INPUT:   COEFFICIENTS  OF  LINEAR  EQUATIONS  FOR  PATH  SEGMENT  AND 
C         REFERENCE  LINE  TO  BE  EVALUATED 

C  OUTPUT:   POINT  OF  INTERSECTION  OF  THE  PATH  SEGMENT  AND  THE  OBSTACLE 
C  REFERENCE  LINE 

IMPLICIT  REAL*8  (A-H,0-Z) 

INTEGER  SEG 

DIMENSION  A1(1000),  B1(1000),  A2(1000),  B2(1000),  D2(1000),  A(2,2) 

A(1,1)  =  A1(LINE) 
A(1,2)  =  B1(LINE) 
A(2,1)  =  A2(SEG) 
A(2,2)  =  B2(SEG) 
BIGEST  =  0.0 

DO  1  L  =  1,2 

DO  1  K  =  1,2 

TEST  =  DABS(A(K,L)) 
IF  (TEST. GT. BIGEST)  THEN 
BIGEST  =  TEST 
KBIG  =  K 
LBIG  =  L 
END  IF 
1  CONTINUE 

IF  (KBIG.EQ.1)  THEN 

IF  (LBIQ.EQ. 1)  THEN 

YINT  =  -D2(SEG)/(B2(SEG)-B1(LINE)*A2(SEG)/A1(LINE)) 

XINT  =  -B1(LINE)*YINT/A1(LINE) 

RETURN 
ELSE 

XINT  =  -D2(SEG)/(A2(SEG)-A1(LINE)*B2(SEG)/B1(LINE)) 

YINT  =  -A1(LINE)*XINT/B1(LINE) 

RETURN 
END  IF 
ELSE 

IF  (LBIG.EQ.1)  THEN 

YINT  =  (D2(SEG)*A1(LINE)/A2(SEG))/ 
+  (B1(LINE)-B2(SEG)*A1(LINE)/A2(SEG)) 

XINT  =  (-D2(SEG)-B2(SEG)*YINT)/A2(SEG) 

RETURN 
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ELSE 

XI  NT  =  (D2(SEG)*B1(LINE)/B2(SEG))/ 
h  (A1(LINE)-A2(SEG)*B1(LINE)/B2(SEG)) 

YINT  =  (-D2(SEG)-XINT*A2(SEG))/B2(SEG) 
RETURN 
END  IF 
END  IF 

END 


SUBROUTINE  ALPHAS  (NEXT,  LOLD,  FSTSTR,  LSTSTR,  SEG,  START,  HEAD, 
+  NSEGS,  PRED) 

C  GIVEN  THAT  A  SEGMENT  OF  THE  PATH  CROSSES  ONLY  ALPHA  RAYS  WHEN 

C  TRAVERSED,  THESE  CROSSINGS  ARE  SIMPLY  INSERTED  INTO  THE  LINK  LIST 

C  ARRAYS  IN  THE  ORDER  IN  WHICH  THEY  WERE  DETECTED,  I.E.,  SMALLEST  TO 

C  LARGEST  IN  ABSOLUTE  VALUE. 

C   INPUT:   STRING  REPRESENTING  THE  REFERENCE  LINES  CROSSED  BY  THE 
C         CURRENT  PATH  SEGMENT,  GIVEN  THAT  ALL  CROSSINGS  ARE  ALPHAS 

C  OUTPUT:   UPDATED  NEXT()  AND  PREDQ  ARRAYS  WHICH  CONTAIN  THE  STRING 
C  REPRESENTING  THE  MOST  CURRENT  PATH  SEGMENT 

INTEGER  NEXT (1000),  FSTSTR,  SEG,  START,  HEAD (1000),  PRED (1000) 

NEXT (LOLD)  =  FSTSTR 
PRED (FSTSTR)  =  LOLD 
LAST  =  LSTSTR  -  1 

IF  ( LAST. GT. FSTSTR)  THEN 
DO  12  I  =  FSTSTR, LAST 
NEXT(I)  =1+1 
PRED(I+1)  =  I 
12     CONTINUE 

ELSE  IF  (LAST. EQ. FSTSTR)  THEN 
NEXT (FSTSTR)  =  LSTSTR 
PRED (LSTSTR)  =  FSTSTR 
END  IF 

IF(SEG.EQ.I)  START  =  LOLD 

LOLD  =  LSTSTR 

RETURN 

END 
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SUBROUTINE  ORDER(SEG,  LINE,  HDINDX,  HEAD,  NEXT,  A1,  A2,  B1,  B2, 
+  D2,  FSTSTR,  LSTSTR,  X,  Y,  START,  LENGTH,  FIRST, 

+  LOLD,  NSEGS,  PRED) 

C  GIVEN  THAT  A  SEGMENT  OF  THE  PATH  CROSSES  ONE  OR  MORE  BETAS,  THIS 

C  SUBROUTINE  DETERMINES  THE  ACTUAL  ORDER  OF  CROSSING  WHEN  TRAVERSING 

C  THE  SEGMENT  FROM  BEGINNING  TO  END.   THIS  IS  DONE  BY  FIRST  DETERMINING 

C  THE  DISTANCE  FROM  THE  SEGMENT  START  POINT  TO  THE  POINT  OF 

C  INTERSECTION  OF  EACH  CROSSED  OBSTACLE  REFERENCE  LINE.   THESE 

C  DISTANCES  ARE  THEN  SORTED  FROM  SMALLEST  TO  LARGEST  AND  CROSSINGS 

C  ARE  UPDATED  IN  THE  LINK  LIST  ACCORDINGLY. 

C   INPUT:   HEAD()  AND  ALL  COORDINATES  REQUIRED  TO  DETERMINE  THE 

C         DISTANCES  FROM  INITIAL  VERTEX  OF  PATH  SEGMENT  TO  EACH  OF 

C         THE  CROSSED  REFERENCE  LINES 

C  OUTPUT:   HEAD(),  NEXT()  AND  PRED()  ARRAYS  CONTAINING  THE  RAW  STRING 
C  WHICH  ACCURATELY  REPRESENTS  THE  PATH  ALONG  THE  CURRENT 

C  SEGMENT 

INTEGER  SEG,  LINE,  HEAD(1000),  NEXT(1000),  FSTSTR,  HDINDX, 
+        START,  FIRST(IOOO),  PRED(1000),  F 

REAL*8  XINTER(IOOO),  YINTER( 1000) ,  A1(1000),  B1(1000), 
+      A2(1000),  B2(1000),  D2(1000),  DIST(1000),  X(1000),  Y(1000) 

DO  8  J  =  FSTSTR,  LSTSTR 
LINE  =  ABS(HEAD(J)) 

CALL  CASES2  (A1,  B1,  A2,  B2,  D2,  SEG,  LINE,  XINTER,  YINTER) 
DIST(J)  =  (XINTER(LINE)-X(SEG))**2+(YINTER(LINE)-Y(SEG))**2 
DIST(J)  =  -DIST(J) 
NEXT(J)  =  J+1 
PRED(J+1)  =  J 
8  CONTINUE 

NEXT (LSTSTR)  =  0 
PRED(LSTSTR+1)  =  LSTSTR 

CALL  MERG2(DIST,  NEXT,  PRED,  FSTSTR,  F) 

NEXT(LOLD)  =  F 
PRED(F)  =  LOLD 
IF(SEG.EQ.I)  START  =  LOLD 
LOLD  =  NEXT (LOLD) 
17  IF(NEXT(LOLD).NE.O)  THEN 

LOLD  =  NEXT (LOLD) 

GO  TO  17 
END  IF 

RETURN 
END 
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SUBROUTINE  CASES2(A1,  B1,  A2,  B2,  D2,  SEG,  LINE,  XINTER,  YINTER) 

C  DETERMINES  COORDINATES  FOR  THE  POINT  OF  INTERSECTION  OF  THE  PATH 

C  SEGMENT  AND  EACH  OF  THE  REFERENCE  LINES  IT  CROSSES.   AGAIN,  IN  ORDER 

C  TO  MAINTAIN  NUMERICAL  STABILITY,  CALCULATIONS  ARE  MADE  USING  THE 

C  MOST  APPROPRIATE  CLOSED  FORM  EQUATION. 

C   INPUT:   COEFFICIENTS  OF  LINEAR  EQUATIONS  FOR  CURRENT  PATH  SEGMENT 
C         AND  ALL  OBSTACLE  REFERENCE  LINES  WHICH  IT  CROSSES 

C  OUTPUT:   COORDINATES  FOR  POINTS  OF  INTERSECTION  OF  EACH  OBSTACLE 
C  REFERENCE  LINE  WITH  THE  PATH  SEGMENT 

IMPLICIT  REAL*8  (A-H,  O-Z) 
INTEGER  SEG 

DIMENSION  A1(1000),  B1(1000),  A2(1000),  B2(1000),  D2(1000), 
+         XINTER(IOOO),  YINTER(IOOO),  A(2,2) 


A(1,1) 

=  A1(LINE) 

A(1,2) 

=  B1(LINE) 

A(2,1) 

=  A2(SEG) 

A(2,2) 

=  B2(SEG) 

BIGEST 

=  0.0 

DO  1  L 

=  1,2 

DO 

1  K  =  1,2 

TEST  =  DABS(A(K,L)) 

IF  (TEST.GT. BIGEST)  THEN 

BIGEST 

=  TEST 

KBIG  = 

K 

LBIG  = 

L 

END  IF 

1  CONTINUE 

IF  (KBIG.EQ.1)  THEN 

IF  (LBIG.EQ.1)  THEN 

YINTER(LINE)=  -D2(SEG)/(B2(SEG)-B1 (L INE)*A2(SEG)/A1 (L INE)) 
XINTER(LINE)=  -B1 (L INE)*YINTER(LINE)/A1(LINE) 
RETURN 
ELSE 

XI NTER(L I NE)=  -D2(SEG)/(A2(SEG)-A1 (L I NE)*B2(SEG)/B1 (L I NE) ) 
YINTER(LINE)=  -A1 (LINE)*XI NTER(LINE)/B1 (LINE) 
RETURN 
END  IF 
ELSE 

IF  (LBIG.EQ.1)  THEN 

YINTER(LINE)  =  (D2(SEG)*A1 (LINE)/A2(SEG))/ 
+  (B1 (L I NE)-B2(SEG)*A1 (L I NE)/A2(SEG) ) 
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XINTER(LINE)  =  (-D2(SEG)-B2(SEG)*YINTER(L INE))/A2(SEG) 
RETURN 
ELSE 

XINTER(LINE)  =  (D2(SEG)*B1 (L INE)/B2(SEG))/ 

(A1(LINE)-A2(SEG)*B1(LINE)/B2(SEG)) 
YINTER(LINE)  =  (-D2(SEG)-XINTER(LINE)*A2(SEG))/B2(SEG) 
RETURN 
END  IF 
END  IF 
END 


SUBROUTINE  MERG2  (DIST,  NEXT,  PRED,  FSTSTR,  F) 

Q  ******************************************************************* 

C  THIS  SUBROUTINE  SORTS  A  SUBSTRING  OF  ALL  POSTIVE  INTEGERS  INTO 
C   INCREASING  ORDER. 

C   INPUT:   A  DOUBLE  LINK  LIST  CONSISTING  OF  DIST,  NEXT,  AND  PRED  ARRAYS 

C  OUTPUT:  A  DOUBLE  LINK  LIST  WITH  ALL  ENTRIES  PLACED  IN  THE  ORDER 
C  IN  WHICH  THEIR  RESPECTIVE  REFERENCE  LINES  WERE  CROSSED 

IMPLICIT  INTEGER(A-Z) 

REAL*8  DIST (1000) 

LOGICAL  DONE 

DIMENSION  NEXT(1000),  PRED(1000) 

DONE  =  .FALSE. 
P  =  FSTSTR 

1  F  =  P 

CALL  SORT2  (F,  PREDF,  P,  DIST,  NEXT,  DONE,  PRED) 

IF (DONE)  RETURN 

GOTO  1 

END 
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SUBROUTINE  SORT2(F,  PREDF,  P,  DIST,  NEXT,  DONE,  PRED) 

C  MARCHES  DOWN  A  SUBSTRING  OF  POSITIVE  INTEGERS  AND  ONLY  SORTS  IF 
C  ELEMENTS  ARE  NOT  IN  INCREASING  ORDER.   IF  A  NUMBER  NEEDS  TO  BE 
C  PLACED  HIGHER  IN  THE  LIST  SUBROUTINE  'PUT'  IS  CALLED  TO  DO  SO 

C   INPUT:   POINTER  F  INTO  DOUBLE  LINK  LIST  ARRAYS  DIST,  NEXT,  AND  PRED 
C         TO  INDICATE  THE  BEGINNING  OF  A  SUBSTING  OF  POSITIVE  INTEGERS, 
C         AND  PREDF 

C  OUTPUT:  P  IS  A  POINTER,  DIST(P)  IS  THE  LAST  POSITIVE  INTEGER  IN  THE 
C         SUBSTRING  THAT  IS  BEGUN  BY  DIST(F) 

IMPLICIT  INTEGER  (A-Z) 
REAL*8  DIST(1000) 
LOGICAL  DONE 

DIMENSION  NEXT(1000),  PRED(1000) 
TAIL  =  F 
1  NTAIL  =  NEXT(TAIL) 

IF(NTAIL.EQ.O)  THEN 

DONE  =  .TRUE. 

RETURN 
END  IF 

IF(DIST(NTAIL).LE.DIST(TAIL))  THEN 

TAIL  =  NTAIL 

GOTO  1 
END  IF 

CALL  PUT2(F,  TAIL,  NTAIL,  PREDF,  DIST,  NEXT,  PRED) 

GOTO  1 

END 
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SUBROUTINE  PUT2  (F,  TAIL,  NTAIL,  PREDF,  DIST,  NEXT,  PRED) 
C  REARRANGES  POINTERS  TO  PLACE  POSITIVE  INTEGERS  IN  INCREASING  ORDER. 


c 

INPUT: 

F 

c 

TAIL 

c 

NTAIL 

c 

c 

c 

PREDF 

c 

OUTPUT: 

SOME  ! 

c 

IN  IT 

START  OF  POSITIVE  SUBSTRING 

END  OF  CURRENTLY  SORTED  PORTION  OF  SUBSTRING 

NEXT(TAIL)  -  POINTER  TO  THE  SUCCESSOR  OF  TAIL  IN  THE 

ITEM  TO  BE  INCORPORATED  INTO  THE  SORTED  PORTION  OF 

THE  LIST 

THE  PREDECESSOR  OF  F  IN  LINKED  LIST 


NTERS  IN  NEXT  AND  PRED  ARE  CHANGED  TO  PUT  DIST (NTAIL) 
TS  PROPER  PLACE  IN  THE  SORTED  PORTION  ON  THE  LIST 


IMPLICIT  INTEGER  (A-Z) 

REAL*8  DIST (1000) 

DIMENSION  NEXT(1000),  PRED(1000) 

IF  (DIST(NTAIL).GE.DIST(F))  THEN 

TEMP  =  NEXT (NTAIL) 

NEXT( NTAIL)  =  F 

PRED(F)  =  NTAIL 

NEXT (TAIL)  =  TEMP 

PRED(TEMP)  =  TAIL 

F  =  NTAIL 

RETURN 
END  IF 

CALL  WALK2(F,  TAIL,  NTAIL,  DIST,  NEXT,  PRED) 
RETURN 
END 
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SUBROUTINE  WALK2(F,  TAIL,  NTAIL,  DIST,  NEXT,  PRED) 

Q  *************************^ 

C  THIS  SUBROUTINE  WALKS  DOWN  THE  LINKED  LIST  AND  PLACES  THE  DIST (NTAIL) 
C   IN  THE  CORRECT  POSITION  IN  THE  DIST  ARRAY. 

C   INPUT:   F,  TAIL,  NTAIL  ARE  USED  AS  ABOVE 

C  OUTPUT:  ALTERS  POINTERS  IN  NEXT  AND  PRED  ARRAYS  TO  PLACE  DIST (NTAIL) 
C         AFTER  DIST(F)  AND  BEFORE  DIST(TAIL)  IN  DOUBLE  LINK  LIST 

IMPLICIT  INTEGER  (A-Z) 
REAL*8  DIST (1000),  HNTAIL 
DIMENSION  NEXT(1000),  PRED(1000) 

I  =  F 

NEXT  I  =  NEXT(I) 

HNTAIL  =  DIST(NTAIL) 

1  IF(HNTAIL.QE.DIST(NEXTI))  THEN 
NEXT(I)  =  NTAIL 
PRED(NTAIL)  =  I 
NNTAIL  =  NEXT(NTAIL) 
NEXT(NTAIL)  =  NEXT  I 
PRED (NEXT  I)  =  NTAIL 
NEXT(TAIL)  =  NNTAIL 
PRED(NNTAIL)  =  TAIL 
RETURN 
END  IF 

I  =  NEXTI 

NEXT!  =  NEXT(I) 

GOTO  1 

END 
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SUBROUTINE  COUNTR( START,  NEXT,  HEAD,  NELEMS) 

Q     *****************************^ 

C  THIS  SUBROUTINE  OOUNTS  THE  NUMBER  OF  ELEMENTS  IN  A  Gl VEN  STR ING  OF 
C  CHARACTERS 

C   INPUT:   HEADQ  AND  NEXT()  ARRAYS  FOR  THE  STRING  OF  CHARACTERS 

C  OUTPUT:   NUMBER  OF  ELEMENTS  IN  THE  STRING 

INTEGER  PTR,  NEXT (1000),  HE AD (1000),  START 

PTR  =  NEXT (START) 
NELEMS  =  0 

20  IF  (PTR.NE.0)  THEN 

NELEMS  =  NELEMS  +  1 
PTR  =  NEXT (PTR) 
GO  TO  20 
END  IF 

RETURN 
END 
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SUBROUTINE  ALG2(NOBS,  NELEMS,  HEAD,  NEXT,  X,  Y,  BX,  BY,  FR,  NR, 
+  PR) 

C  GENERATES  THE  WELL-KNOW  FUNDAMENTAL  GROUP  REPRESENTATION  OF  AN 
C  EQUIVALENCE  CLASS. 

C   INPUT:   RAW  OR  CANONICAL  STRING,  WITH  #  OBSTACLES,  tt  ELEMENTS  IN 
C  STRING,  COORDINATES  OF  A  AND  OBSTACLES  (BK) 

C  OUTPUT:  1)  IF  INPUT  IS  RAW  STRING,  THEN  POSSIBLY  UNREDUCED 
C  FUNDAMENTAL  GROUP  WILL  RESULT  (IF  CANCELLATION  OF  LIKE 

C         POSITIVE  NUMBERS  COULD  HAVE  OCCURRED  IN  THE  RAW  STRING). 
C  2)  IF  INPUT  IS  CANONICAL  STRING,  THEN  THE  RESULTING 

C         FUNDAMENTAL  GROUP  WILL  BE  IN  REDUCED  FORM. 

REAL*8  XA,  YA,  XB,  YB,  X(1000),  Y(1000),  BX(1000),  BY(1000) 
INTEGER  NOBS,  NELEMS,  HEAD(1000),  NEXT(1000),  PRED(1000), 

+        FR(1000),  NR(1000),  PR(1000),  S(1000),  FUNDGP ( 1 000 ) , 

+        START 
LOGICAL  RIGHT(IOOO),  RITE 

LENGTH  =  1 
START  =  NEXT(1) 

XA  =  X(1) 
YA  =  Y(1) 

DO  1  K  =  1,  NOBS 
XB  =  BX(K) 
YB  =  BY(K) 
1    RIGHT(K)  =  RITE(XA,  YA,  XB,  YB) 

NELEM  =  NELEMS  +  1 
DO  10  M  =  2,  NELEM 

S(M)  =  HEAD(START) 

J  =  IABS(S(M)) 

RIGHT(J)  =  .NOT.  RIGHT(J) 

IF(S(M).GT.O)  THEN 

LENGTH  =  LENGTH  +  1 
IF(RIGHT(J))  THEN 

FUNDGP (LENGTH)  =   J 
ELSE 

FUNDGP (LENGTH)  =  -J 
END  IF 
END  IF 

START  =  NEXT (START) 
IF(START.EQ.O)  GOTO  11 
10  CONTINUE 
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11  FR(1)  =  0 
NR(1)  =  2 

DO  2  I  =  2,  LENGTH 

FR(I)  =  FUNDGP(I) 

NR(I)  =  I  +1 

PR(I)  =  I  -  1 
2  CONTINUE 

NR( LENGTH)  =  0 

RETURN 
END 


LOGICAL  FUNCTION  RITE(XA,  YA,  XB,  YB) 

C  THIS  FUNCTION  DETERMINES  WHICH  SIDE  OF  A  GIVEN  DIRECTED  LINE  ANY 
C  POINT  LIES. 

C   INPUT:  TWO  POINTS  THAT  DETERMINE  THE  LINE 

C  OUTPUT:  LOGICAL  VARIABLE  THAT  IS  TRUE  IF  A  POINT  LIES  TO  THE  RIGHT 
C         AND  FALSE  IF  A  POINT  LIES  TO  THE  LEFT 

REAL*8  XB,  YB,  XA,  YA,  SIGMA 

RITE  =  .TRUE. 

SIGNA  =  -((YB*XA)-(XB*YA)) 

IF(SIGNA.LT.O)THEN 

RITE  =  .FALSE. 

RETURN 
END  IF 

IF(SIGNA.EQ.O)  THEN 

PRINT*,'THE  POINT  A  LIES  ON  THE  LINE  LK  PROGRAM  STOPS' 

STOP 
END  IF 

RETURN 
END 
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SUBROUTINE  CANALG2(HEAD,  NEXT,  PRED) 

Q  *************##***:***********************^ 

C  THIS  SUBROUTINE  TAKES  AN  UNREDUCED  FUNDAMENTAL  GROUP  REPRESENTATION 
C  OF  A  GIVEN  CLASS  AND  CANCELS  A  GENERATOR  IF  IT  IS  ADJACENT  TO  ITS 
C   INVERSE. 

C   INPUT:   HEAD,  NEXT,  AND  PRED  ARRAYS 

C  OUTPUT:  HEAD,  NEXT,  AND  PRED  ARRAYS  WITH  NEXT  AND  PRED  REARRANGED 
C         TO  SKIP  AROUND  CANCELLED  ELEMENTS 

IMPLICIT  INTEGER  (A-Z) 

DIMENSION  HEAD(IOOO),  NEXT(1000),  PRED(1000) 

START  =  1 
PTR1  =  START 
PTR2  =  START 
PTR3  =  NEXT(PTR2) 
10  IF  (PTR3.NE.0)  THEN 

IF  (HEAD(PTR2).EQ.-(HEAD(PTR3)))  THEN 
NEXT(PTRI)  =  NEXT(PTR3) 
PRED ( NEXT (PTR3))  =  PTR1 
IF(NEXT(PTR1).EQ.O)  RETURN 
PTR2  =  PTR1 
PTR1  =  PRED(PTRI) 
PTR3  =  NEXT(PTR2) 
ELSE 

PTR1  =  PTR2 
PTR2  =  PTR3 
PTR3  =  NEXT(PTR2) 
END  IF 
GO  TO  10 
END  IF 
RETURN 
END 
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SUBROUTINE  ALG1(HEAD,  NEXT,  PRED,  NELEMS) 

C  THIS  SUBROUTINE  TAKES  A  GIVEN  RAW  STRING  OF  CHARACTERS  REPRESENTING 

C  A  PATH  AND  PRODUCES  THE  CANONICAL  FORM  OF  THAT  STRING.   THIS  IS  DONE 

C  BY  FIRST  ORDERING  ALL  OF  THE  ALPHA  SUBSTRINGS  FROM  SMALLEST  TO 

C  LARGEST  IN  ABSOLUTE  VALUE.   NEXT,  CANCELLATION  IS  PERFORMED  TO 

C  ELIMINATE  ALL  LIKE  PAIRS  OF  ADJACENT  ELEMENTS  FROM  THE  STRING. 

C   INPUT:   RAW  STRING  IN  FORM  OF  DOUBLE  LINKED  LIST  WITH  HEAD(),  NEXTQ 
C         AND  PREDQ  ARRAYS 

C  OUTPUT:   CANONICAL  FORM  OF  THE  RAW  STRING 

IMPLICIT  INTEGER  (A-Z) 

DIMENSION  HEAD(1000),  NEXT(1000),  PRED(1000) 

START  =  1 

CALL  MERG1(HEAD,  NEXT,  PRED) 

CALL  CANCEL (START,  HEAD,  NEXT,  PRED,  NELEMS) 

RETURN 
END 
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SUBROUTINE  MERG1  (HEAD,  NEXT,  PRED) 
C  ************************************************************* 

C  THIS  SUBROUTINE  DOES  THE  INITIAL  SORT  OF  THE  ALPHA  SUBSTRINGS  IN  THE 
C  RAWSTRING  (INCREASING  IN  ABSULUTE  VALUE). 

C   INPUT:   HEAD,  NEXT,  PRED  ARRAYS  REPRESENTING  A  DOUBLE  LINKED  LIST 
C         OF  THE  RAWSTRING  WITH  THE  ALPHA  SUBSTRINGS  UNORDERED 

C  OUTPUT:  POINTERS  STORED  IN  THE  ARRAYS  NEXT  AND  PRED  ARE  ALTERED  SO 

C         THAT  EACH  SUBSTRING  OF  THE  STORED  LIST  WHICH  CONSISTS 

C         ENTIRELY  OF  NEGATIVE  INTEGERS  IS  SORTED  INTO  NON- I NCREAS I NG 

C         ORDER,  WHILE  SUBSTRINGS  OF  POSITIVE  INTEGERS  ARE  LEFT 

C         UNALTERED. 

IMPLICIT  INTEGER(A-Z) 

LOGICAL  DONE 

DIMENSION  HEAD(1000),  NEXT(1000),  PRED(1000) 

DONE  =  .FALSE. 
P  =  1 

1  CALL  FRONT  (P,  F,  PREDF,  HEAD,  NEXT,  DONE) 

IF(DONE)  RETURN 

CALL  SORT  (F,  PREDF,  P,  HEAD,  NEXT,  DONE,  PRED) 

IF (DONE)  RETURN 

GOTO  1 

END 
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SUBROUTINE  FRONT(P,  F,  PREDF,  HEAD,  NEXT,  DONE) 

C  FINDS  THE  BEGINNING  OF  NEGATIVE  INTEGER  STRINGS  (ALPHA  STRING) 

C   INPUT:   POINTER  P  INTO  LINKED  LIST 

C  OUTPUT:  F  IS  A  POINTER  INTO  THE  HEAD  ARRAY.   POINTS  TO  FIRST  NEGATIVE 
C  ENTRY  WHICH  OCCURS  STRICTLY  AFTER  HEAD(P).  PREDF  IS  POINTER 

C         SUCH  THAT  HEAD(PREDF)  IS  THE  PREDESSOR  OF  HEAD(F). 

IMPLICIT  INTEGER(A-Z) 

LOGICAL  DONE 

DIMENSION  HEAD(1000),  NEXT(1000) 

F  =  P 
1  PREDF  =  F 
F  =  NEXT(F) 

IF  (F.EQ.O)  THEN 

DONE  =  .TRUE. 

RETURN 
END  IF 

IF  (HEAD(F).LT.O)  RETURN 

GOTO  1 

END 
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SUBROUTINE  SORT(F,  PREDF,  P,  HEAD,  NEXT,  DONE,  PRED) 

C  MARCHES  DOWN  A  SUBSTRING  OF  NEGATIVE  INTEGERS  AND  ONLY  SORTS  IF 
C  ELEMENTS  ARE  IN  INCREASING  ORDER.   IF  A  NUMBER  NEEDS  TO  BE 
C  PLACED  HIGHER  IN  THE  LIST  SUBROUTINE  'PUT'  IS  CALLED  TO  DO  SO 

C   INPUT:   POINTER  F  INTO  DOUBLE  LINK  LIST  ARRAYS  HEAD,  NEXT,  AND  PRED 
C         TO  INDICATE  THE  BEGINNING  OF  A  SUBSTRING  OF  NEGATIVE 
C         INTEGERS;  AND  PREDF 

C  OUTPUT:  PISA  POINTER,  HEAD(P)  IS  THE  LAST  NEGATIVE  INTEGER  IN  THE 
C  SUBSTRING  THAT  IS  BEGUN  BY  HEAD(F) 

IMPLICIT  INTEGER  (A-Z) 

LOGICAL  DONE 

DIMENSION  HEAD(1000),  NEXT(1000),  PRED(1000) 

TAIL  =  F 
1  NTAIL  =  NEXT(TAIL) 

IF(NTAIL.EQ.O)  THEN 

DONE  =  .TRUE. 

RETURN 
END  IF 

IF(HEAD(NTAIL).GT.O)  THEN 

P  =  TAIL 

RETURN 
END  IF 

IF(HEAD(NTAIL).LE.HEAD(TAIL))  THEN 

TAIL  =  NTAIL 

GOTO  1 
END  IF 

CALL  PUT(F,  TAIL,  NTAIL,  PREDF,  HEAD,  NEXT,  PRED) 

GOTO  1 

END 
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SUBROUTINE  PUT  (F,  TAIL,  NTAIL,  PREDF,  HEAD,  NEXT,  PRED) 
C  ************************************************************ 

C  SUBROUTINE  THAT  REARRANGES  POINTERS  TO  PLACE  A  NEGATIVE  INTEGER  IN 
C  NON- INCREASING  ORDER. 

C  INPUT:   F     -  START  OF  NEGATIVE  SUBSTRING 

C         TAIL-END  OF  CURRENTLY  SORTED  PORTION  OF  SUBSTRING 

C  NTAIL  =  NEXT(TAIL)  -  POINTER  TO  THE  SUCCESSOR  OF  TAIL  IN  THE 

C  ITEM  TO  BE  INCORPORATED  INTO  THE  SORTED  PORTION  OF 

C  THE  LIST. 

C         PREDF  -  THE  PREDECESSOR  OF  F  IN  LINKED  LIST 

C  OUTPUT:  SOME  POINTERS  IN  NEXT  AND  PRED  ARE  CHANGED  TO  PUT  HEAD(NTAIL) 
C         IN  ITS  PROPER  PLACE  IN  THE  SORTED  PORTION  ON  THE  LIST 

IMPLICIT  INTEGER  (A-Z) 

DIMENSION  HEAD(1000),  NEXT(1000),  PRED(1000) 

IF  (HEAD(NTAIL).GE.HEAD(F))  THEN 

NEXT (PREDF)  =  NTAIL 

PRED (NTAIL)  =  PREDF 

TEMP  =  NEXT (NTAIL) 

NEXT(NTAIL)  =  F 

PRED(F)  =  NTAIL 

NEXT (TAIL)  =  TEMP 

PRED (TEMP)  =  TAIL 

F  =  NTAIL 

RETURN 
END  IF 

CALL  WALK(F,  TAIL,  NTAIL,  HEAD,  NEXT,  PRED) 

RETURN 
END 
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SUBROUTINE  WALK(F,  TAIL,  NTAIL,  HEAD,  NEXT,  PRED) 

C  THIS  SUBROUITNE  WALKS  DOWN  THE  LINKED  LIST  AND  PLACES  THE  HEAD (NTAIL) 
C   IN  THE  CORRECT  POSITION  IN  THE  HEAD  ARRAY.  (DECREASING  ORDER) 

C   INPUT:   F,  TAIL,  NTAIL  ARE  USED  AS  ABOVE 

C  OUTPUT:  ALTERS  POINTERS  IN  NEXT  AND  PRED  ARRAYS  TO  PLACE  HEAD (NTAIL) 
C         AFTER  HEAD(F)  AND  BEFORE  HEAD (TAIL)  IN  DOUBLE  LINK  LIST 

IMPLICIT  INTEGER  (A-Z) 

DIMENSION  HEAD(1000),  NEXT(1000),  PRED(1000) 

I  =  F 

NEXT  I  =  NEXT(I) 

HNTAIL  =  HEAD(NTAIL) 

1  IF(HNTAIL.G£.HEAD(NEXTI))  THEN 
NEXT(I)  =  NTAIL 
PRED( NTAIL)  =  I 
NNTAIL  =  NEXT( NTAIL) 
NEXT (NTAIL)  =  NEXT  I 
PRED (NEXT  I)  =  NTAIL 
NEXT(TAIL)  =  NNTAIL 
PRED(NNTAIL)  =  TAIL 
RETURN 
END  IF 

I  =  NEXT  I 

NEXT  I  =  NEXT(I) 

GOTO  1 

END 


83 


SUBROUTINE  CANCEL (START,  HEAD,  NEXT,  PRED,  NELEMS) 
q   **************************^ 

C  GIVEN  THE  SORTED  RAW  STRING,  THIS  SUBROUTINE  REDUCES  THE  STRING  TO 

C  CANONICAL  FORM  BY  ADJUSTING  THE  LINK  LIST  TO  SKIP  ANY  PAIRS  OF 

C  ADJACENT  LIKE  ELEMENTS  IN  THE  STRING.   ONCE  TWO  ELEMENTS  ARE 

C  ELIMINATED  FROM  THE  STRING,  THE  SUBSTRINGS  WHICH  WERE  ON  THEIR  LEFT 

C  AND  RIGHT  CONCATENATE  TO  FORM  A  NEW  STRING  AND  HENCE  A  NEW  PAIR  OF 

C  ADJACENT  ELEMENTS  WHICH  MUST  ALSO  BE  CHECKED  FOR  EQUALITY.   IN  THE 

C  EVENT  A  BETA  SUBSTRING  BECOMES  ANNHILATED  FROM  THE  STRING,  THE  TWO 

C  ADJACENT  ALPHA  STRINGS  CONCATENATE  AND  ARE  AGAIN  SORTED  FROM  SMALLEST 

C  TO  LARGEST  IN  ABSOLUTE  VALUE  AS  A  SINGLE  SUBSTRING. 

C   INPUT:   SORTED  RAW  STRING 

C  OUTPUT:   CANONICAL  FORM  OF  THE  RAW  STRING 

IMPLICIT  INTEGER  (A-Z) 

DIMENSION  HEAD(1000),  NEXT(1000),  PRED(1000) 

START  =  1 
START 1  =  START 
START2  =  START 
BETA  =  START 
PTR1  =  START 1 
PTR2  =  START 1 
50  PTR3  =  NEXT(PTR2) 

IF  (HEAD(PTR3).LT.O)  THEN 

START 1  =  PTR3 

START2  =  START 1 
ELSE 

BETA  =  PTR3 
END  IF 

10  IF  (PTR3.NE.0)  THEN 

CALL  CHEKER  (HEAD,  START1,  START2,  PTR1,  PTR2,  PTR3,  BETA) 
IF  (HEAD(PTR2).EQ.HEAD(PTR3))  THEN 
NEXT(PTRI)  =  NEXT(PTR3) 
PRED(NEXT(PTR3))  =  PTR1 
PTR2  =  NEXT(PTR3) 
IF  (PTR2.EQ.0)  GO  TO  60 
CHECK  1  =  HEAD(PTR2)  *  HEAD(PTR3) 
CHECK2  =  HEAD(PTRI)  *  HEAD(PTR2) 

IF  (CHECK 1.LT.0. AND. CHECK2.GT.0. AND. HEAD (PTR2).LT.0)THEN 
PTR2  =  PTR1 
PTR1  =  PRED(PTRI) 
PTR3  =  NEXT(PTR2) 
START 1  =  START2 
6  IF(START1.NE.1.AND.HEAD(PRED(START1)).LT.0)THEN 
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START 1  =  PRED ( START 1) 
GOTO  6 
END  IF 

START2  =  PTR3 

CALL  MERGE  (HEAD,  NEXT,  START 1,  START2,  PRED) 
PTR1  =  PRED ( START 1) 
PTR2  =  PTR1 
QO  TO  50 
ELSE  IF  (CHECK1.LT.O.AND.CHECK2.EQ.O.AND.HEAD(PTR2).LT.O) 
+        THEN 

PTR2  =  PTR1 
PTR3  =  NEXT(PTRI) 
START 1  =  PTR3 
START2  =  START 1 
00  TO  10 
ELSE  IF  (CHECK1.GT.O.AND.CHECK2.GT.O.AND.HEAD(PTR2).GT.O) 
+        THEN 

IF  (BETA.EQ.2)  THEN 

BETA  =  1 
ELSE  IF  (BETA.EQ.1)  THEN 

00  TO  9 
ELSE 

IF(STARTI.NE.I)  BETA  =  PRED(STARTI) 
11  IF(HEAD(PRED(BETA)).GT.O)  THEN 

BETA  =  PRED (BETA) 
GO  TO  11 
END  IF 
END  IF 
9  PTR1  =  PRED (BETA) 

PTR2  =  BETA 
PTR3  =  NEXT(PTR2) 
GO  TO  10 
ELSE  IF  (CHECK1.LT.O.AND.CHECK2.GT.O.AND.HEAD(PTR2).GT.O) 
+        THEN 

PTR2  =  PTR1 

PTR1  =  PRED(PTRI) 

IF(HEAD(PTR1).LT.0)THEN 

START 1  =  PTR1 

5  IF  ( HEAD ( PRED ( START 1)).LT.0)THEN 

START 1  =  PRED ( START 1) 
GO  TO  5 
END  IF 
ENDIF 

PTR3  =  NEXT(PTR2) 
GO  TO  10 
ELSE  IF  (CHECK1.LT.O.AND.CHECK2.LT.O.AND.HEAD(PTR2).LT.O) 
+        THEN 

IF(START1.EQ.1)  THEN 
START 1  =  PTR2 
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START2  =  PTR2 
PTR3  =  NEXT(PTR2) 
ELSE 

START 1  =  START2 
PTR3  =  NEXT(PTR2) 
END  IF 
GO  TO  10 
ELSE 

START 1  =  START2 
PTR3  =  NEXT(PTR2) 
END  IF 
ELSE 

PTR1  =  PTR2 
PTR2  =  PTR3 
PTR3  =  NEXT(PTR2) 
END  IF 
00  TO  10 
END  IF 

60  CALL  COUNTR  (START,  NEXT,  HEAD,  NELEMS) 

RETURN 
END 
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SUBROUTINE  CHEKER(HEAD,  START 1 ,  START2,  PTR1,  PTR2,  PTR3,  BETA) 

Q  ***********************^^ 

C  AS  POINTERS  MOVE  ALONG  THE  CHARACTER  STRING  DURING  CANCELLATION  THIS 

C  SUBROUTINE  CHECKS  TO  DETERMINE  WHETHER  OR  NOT  THE  END  OF  ONE 

C  SUBSTRING  IS  REACHED  AND  A  NEW  ONE  BEGINS.   DEPENDING  ON  THE  OUTCOME 

C  OF  THIS  CHECK,  THE  POINTERS  USED  TO  IDENTIFY  THE  BEGINNING  OF  THE 

C  TWO  LATEST  ALPHA  STRINGS  AND  THE  LATEST  BETA  STRING  ARE  UPDATED 

C   INPUT:   HEAD()  ARRAY  AND  PRESENT  POINTER  LOCATIONS  FROM  'CANCEL' 
C         SUBROUT I NE 

C  OUTPUT:   MODIFIED  INDICES  FOR  LOCATION  OF  ALPHA  AND  BETA  STRINGS 

INTEGER  HEAD(1000),  START 1 ,  START2,  PTR1,  PTR2,  PTR3,  CHECK,  BETA 

CHECK  =  HEAD(PTR3)  *  HEAD(PTR2) 

IF  (CHECK. GT.O)  THEN 

IF  (HEAD(PTR2).LT.O. AND. HEAD(PTRI). GT.O)  START2  =  PTR2 
IF  (HEAD(PTR2). GT.O. AND. HEAD(PTRI).LT.O)  BETA  =  PTR2 
ELSE  IF  (CHECK. LE.O)  THEN 

IF  (HEAD(PTR2). GT.O. AND. HEAD(PTRI).LT.O)  THEN 
BETA  =  PTR2 
START 1  =  START2 
START2  =  PTR3 
ELSE  IF  (HEAD(PTR2).GT.O)  THEN 
IF  (START1.NE.1)  THEN 
START 1  =  START2 
START2  =  PTR3 
ELSE 

START 1  =  PTR3 
START2  =  START 1 
END  IF 
ELSE  IF  (HEAD(PTR2).LT.O. AND. HEAD(PTRI). GT.O)  THEN 
BETA  =  PTR3 
START 1  =  START2 
START2  =  PTR2 
ELSE  IF  (HEAD(PTR2).LT.O.AND.HEAD(PTR1).LT.O)  THEN 

BETA  =  PTR3 
END  IF 
END  IF 

RETURN 
END 
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SUBROUTINE  MERGE (HEAD,  NEXT,  START 1 ,  START2,  PRED) 

C  THIS  SUBROUTINE  MERGES  TWO  ORDERED  ALPHA  STRINGS.   THE  RESULTING 
C  SUBSTRING  WILL  BE  INCREASING  IN  ABSOLUTE  VALUE. 

C  INPUT:  HEADQ,  NEXT()  AND  PRED()  ARRAYS  ALONG  WITH  INDICES  START1 
C  AND  START2  INDICATING  BEGINNING  OF  THE  ALPHA  STRINGS  TO  BE 
C         MERGED 

C  OUTPUT:  MODIFIED  STRING  CONTAINING  THE  MERGED  ALPHA  STRINGS  WHICH 
C  WERE  A  RESULT  OF  ANNHILATION  OF  A  BETA  STRING 

INTEGER  HEAD(1000),  NEXT(1000),  START 1 ,  START2,  P1,  P2,  TAIL, 
+  PRED(1000) 

LOGICAL  GO 

IF  (START1.EQ.1)  GO  TO  2 

IF(HEAD(PRED(START1)).EQ.O)  START1  =  NEXT ( PRED (START1 )) 

IF(NEXT(PRED(START1)).NE.START1)THEN 
START1  =  NEXT(PRED(START1)) 

7  IF  (START1.EQ.1)  GO  TO  2 

I F(HEAD(START1 ) .GT.0)THEN 

START 1  =  NEXT ( START 1) 

GO  TO  7 
END  IF 

8  IF(HEAD(PRED(START1)).LT.0)THEN 

START 1  =  PRED ( START 1) 
GO  TO  8 
END  IF 
END  IF 

2  CALL  SET(HEAD,  NEXT,  START1,  START2,  P1,  P2,  TAIL,  PRED) 

1  IF(GO(P1,  P2,  START2,  HEAD,  NEXT))  THEN 

CALL  MERGER (HEAD,  NEXT,  P1,  P2,  TAIL,  PRED) 
GOTO  1 
END  IF 

CALL  FXTAIL(HEAD,  NEXT,  P1,  P2,  TAIL,  START2,  PRED) 

RETURN 
END 


88 


SUBROUTINE  SET (HEAD,  NEXT,  START 1 ,  START2,  P1,  P2,  TAIL,  PRED) 

Q  ******************************** 

C  SETS  THE  STARTING  POINTER  OF  THE  NEW  STRING  AND  THE  TWO  POINTERS 
C  OF  THE  STRINGS  TO  BE  MERGED 

C   INPUT:   POINTERS  INTO  TWO  ALPHA  STRINGS 

C  OUTPUT:  POINTER  INTO  BEGINNING  OF  A  NEW  SORTED  STRING 

INTEGER  HEAD(1000),  NEXT(1000),  START 1 ,  START2,  P1,  P2,TAIL, 
+        PRED (1000) 

P1  =  START 1 
P2  =  START2 
IF(HEAD(P1).GT.HEAD(P2))  THEN 

TAIL  =  P1 

START 1  =  P1 

P1  =  NEXT(P1) 
ELSE 

TAIL  =  P2 

START 1  =  P2 

NEXT(PRED(P1))  =  P2 

PRED(P2)  =  PRED(P1) 

P2  =  NEXT(P2) 
END  IF 

RETURN 
END 
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LOGICAL  FUNCTION  QO(P1,  P2,  START2,  HEAD,  NEXT) 

C  LOGICAL  FUNCTION  TO  DETERMINE  WHEN  THE  END  OF  EITHER  LIST  IS  REACHED 

C   INPUT:   POINTERS  INTO  TWO  ALPHA  SUBSTRINGS 

C  OUTPUT:  LOGICAL  VARIABLE  WHICH  IS  TRUE  IF  EITHER  POINTER  IS  AT  THE 
C  END  OF  A  STRING,  AND  FALSE  OTHERWISE 

INTEGER  HEAD(1000),  NEXT(1000),  P1,  P2,  START2 

GO  =  .TRUE. 

IF(P1.EQ.START2)THEN 

GO  =  .FALSE. 

RETURN 
END  IF 

IF(NEXT(P2).EQ.0)THEN 

GO  =  .FALSE. 

RETURN 
END  IF 

IF(HEAD(P2).GT.0)THEN 

GO  =  .FALSE. 

RETURN 
END  IF 

RETURN 
END 
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SUBROUTINE  MERGER(HEAD,  NEXT,  P1,  P2,  TAIL,  PRED) 

C  THIS  SUBROUTINE  COMPARES  TWO  ELEMENTS  -  ONE  IN  EACH  ALPHA  SUBSTRING 
C  AND  PLACES  THE  ONE  THAT  IS  THE  SMALLEST  ON  THE  LIST  REPRESENTING  THE 
C  NEW  SORTED  LIST 

C   INPUT:   HEAD,  NEXT,  PRED,  AND  POINTERS  INTO  THE  HEAD  ARRAY 

C  OUTPUT:  HEAD,  NEXT,  PRED,  AND  POINTERS  FOR  THE  NEXT  TWO  ELEMENTS 
C         TO  BE  CHECKED 

INTEGER  HEAD(1000),  NEXT(1000),  P1,  P2,  TAIL  ,  PRED(1000) 

IF(HEAD(P1).GT.HEAD(P2))  THEN 

CALL  TACON1(NEXT,  P1,  TAIL,  PRED) 
ELSE 

CALL  TACON2(NEXT,  P2,  TAIL,  PRED) 
END  IF 

RETURN 
END 


SUBROUTINE  TACON1  (NEXT,  P1,  TAIL,  PRED) 

C  ADDS  THE  ELEMENT  TO  WHICH  P1  POINTS  TO  THE  TAIL  OF  THE  NEW  STRING 

C   INPUT:   NEXT,  PRED,  P1,  TAIL 

C  OUTPUT:  NEXT,  PRED,  P1,  TAIL 

INTEGER  NEXT(1000),  P1,  TAIL,  PRED(1000) 

NEXT(TAIL)  =  P1 
PRED(P1)  =  TAIL 
TAIL  =  P1 
P1  =  NEXT(P1) 

RETURN 
END 
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SUBROUTINE  TAOON2  (NEXT,  P2,  TAIL,  PRED) 

C  ADDS  THE  ELEMENT  TO  WHICH  P2  POINTS  TO  THE  TAIL  OF  THE  NEW  STRING 

C   INPUT:   NEXT,  PRED,  P2,  TAIL 

C  OUTPUT:  NEXT,  PRED,  P2,  TAIL 

INTEGER  NEXT(1000),  P2,  TAIL,  PRED(1000) 

NEXT(TAIL)  =  P2 

PRED(P2)  =  TAIL 

TAIL  =  P2 

P2  =  NEXT(P2) 

RETURN 

END 
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SUBROUTINE  FXTAI L(HEAD,  NEXT,  P1,  P2,  TAIL,  START2,  PRED) 

C  ATTACHES  THE  END  OF  THE  LONGER  SORTED  STRING  TO  THE  END  OF  THE  MERGED 
C  STR I NG 

C   INPUT:   P1,  P2,  HEAD,  NEXT,  PRED,  TAIL 

C  OUTPUT:  P1,  P2,  HEAD,  NEXT,  PRED,  TAIL,  START2 

INTEGER  HEAD(1000),  NEXT(1000),  PI,  P2,  TAIL,  START2,  P11, 
+        PRED(1000) 

IF(P1.EQ.START2)THEN 

CALL  TACON2(NEXT,  P2,  TAIL,  PRED) 

RETURN 
END  IF 

IF(NEXT(P2).EQ.0)THEN 

CALL  TACON1(NEXT,  P1,  TAIL,  PRED) 
P11  =  TAIL 

2  IF(P1.NE.START2)THEN 

P11  =  P1 
P1  =  NEXT(P1) 
GOTO  2 
END  IF 

NEXT(P11)  =  P2 
PRED(P2)  =  P11 
RETURN 
END  IF 

IF(HEAD(P2).GT.0)THEN 

CALL  TACON1(NEXT,  P1,  TAIL,  PRED) 
P11  =  TAIL 

3  IF(P1.NE.START2)  THEN 

P11  =  P1 

P1  =  NEXT(P1) 

GOTO  3 
END  IF 

NEXT(P11)  =  P2 
PRED(P2)  =  P11 
RETURN 
END  IF 

RETURN 
END 
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SUBROUTINE  CTREX(FL,  NL,  PL,  FR,  NR,  PR,  BSEED,  PSEED) 
C  ********************************************************** 

C  THIS  SUBROUTINE  CHECKS  F(C(R(P)))  AGAINST  K(F(R(P)))  FOR  EQUALITY. 

C   INPUT:   FUNDAMENTAL  GROUP  REPRESENTATIONS  FOR  BOTH  ALGORITHMS 

C  OUTPUT:   MESSAGE  INDICATING  THE  OCCURRENCE  OF  UNEQUAL  FUNDAMENTAL 
C  GROUPS 

INTEGER  FL(1000),  NL(1000),  PL(1000),  FR(1000),  NR(1000),  PR(1000) 

IF(NL(1).EQ.O.AND.NR(1).EQ.O)  GOTO  11 

I  =  NL(1) 
J  =  NR(1) 

10  IF  ( I .EQ.0.AND.J.EQ.0)  GO  TO  11 

IF  (FL(l).EQ.O.AND.FR(J).EQ.O)  GO  TO  11 

IF  (FL(I).EQ.FR(J))  THEN 
I  =  NL(I) 
J  =  NR(J) 
GOTO  10 
ELSE 

PR  I  NT* , ' STR I NGS  UNEQUAL ' 
PR  I  NT* , '  ' 
IF(NL(1).EQ.0)THEN 

PR  I  NT*,' THE  GROUP  IS  EMPTY  AFTER  ALG1  SEQUENCE' 
ELSE 

PR  I  NT*,' THE  FUNDAMENTAL  GROUP  AFTER  ALG1  IS:' 

CALL  PRINTS(FL,NL) 
END  IF 
IF(NR(1).EQ.O)  THEN 

PR  I  NT* , ' THE  GROUP  I S  EMPTY  AFTER  ALG2  SEQUENCE ' 
ELSE 

PRINT*, 'THE  FUNDAMENTAL  GROUP  AFTER  ALG2  IS:' 

CALL  PRINTS(FR,NR) 
END  IF 
PRINT*,'  ' 
PR  I  NT* ,' BSEED  AND  PSEED  ARE' , BSEED, PSEED 

END  IF 

11  RETURN 
END 
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SUBROUTINE  PR  I  NTS (HEAD,  NEXT) 
C  ******************************************************** 

C  THIS  SUBROUTINE  CAN  BE  USED  TO  PRINT  OUT  ANY  STRING  THAT  IS  STORED  IN 
C  LINK  LIST  FORM 

C   INPUT:   LINK  LIST  ARRAYS 

C  OUTPUT:  HORIZONTAL  STRING  OF  THE  ELEMENTS  ACCORDING  TO  THE  NEXT() 

INTEGER  HEAD(1000),  NEXT(1000),  STRING(IOOO) 

NSTART  =  NEXT(1) 
I  =  0 

98  IF(NSTART.NE.O)  THEN 
1  =  1  +  1 

STRING(I)  =  HEAD(NSTART) 
NSTART  =  NEXT (NSTART) 
GOTO  98 
END  IF 

PRINT  111,  (STRING(J),  J  =  1,  I) 
111  FORMATC  ',2014) 
PRINT*,'  ' 

RETURN 
END 
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C  ALL  SUBROUTINES  BELOW  THIS  POINT  WERE  USED  AS  DEBUG  TOOLS.   THEY 

C  GENERATE  CRUDE  GRAPHS  WHICH  PLOT  THE  BOARDS  AND  THE  PATHS. 

C  THE  VALUES  AQUIRED  BY  THE  RAWSTRING  ,  MERGE,  AND  MANY  OTHER 

C  SUBROUTINES  WERE  CHECKED  USING  THESE  ROUTINES. 

SUBROUTINE  GRAPH (BX,  BY,  X,  Y,  NOBS,  NUMPTS) 

REAL*8  X(1000),  Y(1000),  BX(1000),  BY(1000) 
CHARACTER*!  MATR(53,  105) 

CALL  GRAFPA(X,  Y,  NUMPTS,  MATR) 
CALL  GRAFOB(BX,  BY,  NOBS,  MATR) 

DO  1  I  =  1,53 
1    PRINT  111,  (MATR(I,J),  J=1,105) 
111   FORMATC  *,105A1) 

RETURN 
END 
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SUBROUTINE  OR AFP A   (X,    Y,    NUMPTS,    MATR) 

C.  ^r^h^h  ^p  ^p  ^p^n  ^p  ^p  ^p  ^p  ^p  ^*  ^P  ^P  ^P  ^p  ^p  *T*  ^p  ^p  ^p  ^p  ^*  ^p  *^^p^p  ^p^p  ^p  ^p  ^p  ^p  ^p^»^p^p^p^p*^^p^p^p  ^p  ^p  ^p^p  «^^p^p^^^^^p^p^^p^P^^^^^P^P^P^P^P^^^P^P^P^P^^ 

REAL*8  X(1000),    Y(1000) 
CHARACTER*1    MATR(53,    105),    CH 

DO   10  L   =   2,52 

DO   10  M  =   2,104 
10  MATR(L,M)    =    '     ' 

DO   1    I    =    1,    105 

MATR(1,I)    =    **' 

1  MATR(53, I)   =    '*' 

DO  4    I    =    1 ,    53 

MATRO  ,1)    =    '*' 
4  MATR(I ,105)   =    '*' 

XK   =   X(1) 

YK  =  Y(1) 

CALL  COOR(XK,    YK,     IX,    IY) 

CALL  CHARPA(1,    CH) 

MATR(IY, IX)    =  CH 

DO  2  K   =   2, NUMPTS 

2  CALL   FILL(K,    X,    Y,    MATR) 

DO  3  K   =   2, NUMPTS 
XK   =   X(K) 
YK   =   Y(K) 

CALL  COOR(XK,    YK,     IX,     IY) 
CALL  CHARPA(K,    CH) 

3  MATR(IY, IX)   =  CH 

RETURN 
END 


SUBROUTINE  COOR(BX,    BY,     IX,     IY) 

REAL*8   BX,    BY 

IX   =  MESHP(BX)   *   2 
IY   =   54   -   MESHP(BY) 

RETURN 
END 
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INTEGER  FUNCTION  MESHP(X) 

REAL*8  X 

MESHP  =  ((X+1.)*(51./2.))  +  2. 

RETURN 
END 


SUBROUTINE  CHARPA(K,  CH) 

CHARACTER* 1  CH 

IF(K.EQ.I)  THEN 

CH  =  'A' 

RETURN 
END  IF 

IF(K.EQ.2)  THEN 

CH  =  'B' 

RETURN 
END  IF 

IF(K.EQ.3)  THEN 

CH  =  'C 

RETURN 
END  IF 

IF(K.EQ.4)  THEN 

CH  =  'D' 

RETURN 
END  IF 

IF(K.EQ.5)  THEN 

CH  =  'E* 

RETURN 
END  IF 

IF(K.EQ.6)  THEN 

CH  =  'A' 

RETURN 
END  IF 

END 
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SUBROUTINE  FILL(K,  X,  Y,  MATR) 

REAL*8  X(1000),  Y(1000),  U,  V,  DX,  DY,  XK,  YK,  XKM1 ,  YKM1 
CHARACTER*1  MATR(53,105) 

XK  =  X(K) 
YK  =  Y(K) 
XKM1  =  X(K-1) 
YKM1  =  Y(K-1) 

CALL  O0OR(XK,  YK,  KX,  KY) 

CALL  OOOR(XKM1,  YKM1 ,  KXM1 ,  KYM1) 

L  =  IABS(KX  -  KXM1)  -  1 
M  =  IABS(KY  -  KYM1)  -  1 
IF  (M.GT.L)  L  =  M 
DX  =  (XK  -  XKM1)/(L+1) 
DY  =  (YK  -  YKM1)/(L+1) 
DO  1  J  =  1,L 

U  =  XKM1  +  J*DX 

V  =  YKM1  +  J*DY 

CALL  COOR(U,V, IU, IV) 
1     MATR(IV,  IU)  =  '+' 

RETURN 
END 


SUBROUTINE  GRAFOB  (X,  Y,  NOBS,  MATR) 

REAL*8  X(1000),  Y(1000) 
CHARACTER* 1  MATR (53,  105),  CH 

MATR(27,  53)  =  '*' 

DO  2  K  =  1,NOBS 
XK  =  X(K) 
YK  =  Y(K) 

CALL  COOR(XK,  YK,  IX,  IY) 
CALL  CHAR(K,  CH) 
2     MATR(IY,IX)  =  CH 

RETURN 
END 
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SUBROUTINE  CHAR(K,  CH) 
CHARACTER* 1  CH 

IF(K.EQ.1)  THEN 

CH  =  '1' 

RETURN 
END  IF 

IF(K.EQ.2)  THEN 

CH  =  '2' 

RETURN 
END  IF 

IF(K.EQ.3)  THEN 

CH  =  '3' 

RETURN 
END  IF 

IF(K.EQ.4)  THEN 

CH  =  '4' 

RETURN 
END  IF 

IF(K.EQ.5)  THEN 
CH  =  '5' 

RETURN 
END  IF 

IF(K.EQ.6)  THEN 

CH  =  *6' 

RETURN 
END  IF 

IF(K.EQ.7)  THEN 

CH  =  '7' 

RETURN 
END  IF 

IF(K.EQ.8)  THEN 

CH  =  '8' 

RETURN 
END  IF 

IF(K.EQ.9)  THEN 

CH  =  '9' 

RETURN 
END  IF 

IF(K.EQ.IO)  THEN 
CH  =  '©' 
RETURN 
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END  IF 

IF(K.EQ.11)  THEN 

CH  =  '#' 

RETURN 
END  IF 

IF(K.EQ.12)  THEN 

CH  =  '$' 

RETURN 
END  IF 

IF(K.EQ.13)  THEN 

CH  =  '%' 

RETURN 
END  IF 

IF(K.EQ.U)  THEN 

CH  =  '&' 

RETURN 
END  IF 

IF(K.EQ.15)  THEN 

CH  =  V 

RETURN 
END  IF 

IF(K.EQ.16)  THEN 

CH  =  '/' 

RETURN 
END  IF 

IF(K.EQ.17)  THEN 

CH  =  'V 

RETURN 
END  IF 

IF(K.EQ.18)  THEN 

CH  =  '<' 

RETURN 
END  IF 

IF(K.EQ.19)  THEN 

CH  =  '>' 

RETURN 
END  IF 

IF(K.EQ.20)  THEN 

CH  =  '?' 

RETURN 
END  IF 
END 
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APPENDIX  B.   A  SMALL  EXAMPLE 

Given  the  following  five  classes  (pictured  in  Figure  B.l) 
as  candidates  for  the  class  containing  the  true  shortest  path: 
We  formulate  a  lower  bound  for  each  class. 

1  .      OtiCt* 

2  .  ajB^ 

3  .  ct,B2 

4  .  ajB^o^Bj 
5.  a^ctiBi 


Figure  B.l   The  Five  Classes  Associated  With  Two  Obstacles 
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Consider  the  class  atifi^  pictured  below  in  Figure  B.2 


a 


2  "^  "1 

Figure  B.2   The  Class  a^ 


To  formulate  the  lower  bound  for  each  class  we  first 
shrink  the  obstacles  down  to  a  point  (Figure  B.3).  In  order 
to  cross  ctj  we  develop  a  cone  of  directions  (Figure  B.3)  in 
which  the  first  segment  of  the  bounding  polygonal  path  must 
lie.  „  P9 


a 


Figure  B.3   The  Cone  of  Directions  Associated  With  a 
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This  cone  is  then  intersected  with  the  one  associated  with  an 
Cj  crossing  pictured  in  Figure  B.4. 


Figure  B.4   Cone  Associated  With  B 


The  resulting  cone  of  directions  is  denoted  by  the  double  arc 
in  Figure  B . 5 . 

02 


OCi 


Figure  B.5   The  Intersected  Cone 
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Since  we  have  reached  the  end  of  the  class  name  we  check 
to  see  if  z  lies  in  this  cone  of  directions.  If  z  lies  in  the 
cone  we  go  directly  to  z.  Otherwise,  we  determine  which  side 
of  the  cone  z  lies  so  that  we  can  find  the  peg  on  which  the 
path  bends . 

In  this  example,  z  does  not  lie  in  the  final  cone  of 


directions  so  the  path  bends  at  the  point  b 


Thus  the 


bounding  path  is  made  of  two  line  segments.   One  from  a  to  bj 
and  the  other  from  bj  to  z  (Figure  B.6). 


a 


a 


-2'  -  -1 

Figure  B.6   The  Bounding  Path  for  the  Class  a^ 


A  similar  procedure  is  done  on  each  of  the  remaining  four 
classes  yielding  the  following  lower  bounds: 

1.5231 
1.5303 
2.2885 
3.2885 
2.7008 


1 .  a^a 

2  .  a2fi 

3  .  Qjfl 

4  .  ^Q 


5  .    ctjfi 


Otjftj 


aft 
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This  list  is  then  put  into  increasing  order: 

1.5231 

1.5303 

2.2885 

2.7008 

3.2885 

Now  the  first  class  on  this  ordered  list  is  chosen.  To 
find  the  true  shortest  path  we  expand  the  obstacles  back  to 
their  original  shape.  If  the  path  bends  on  a  peg,  then  a 
tangent  line  is  constructed  from  a  start  point  to  the  correct 
side  of  the  obstacle  (Figure  B.7). 


1 .  ct^ 

2  .  a22>l 

3  .  0:^2 

4.  a^a^Bj 

5  .  ctjBi  aoB2 


a2^  ^  ^ 

Figure  B.7   A  Tangent  Line  From  a  to  B2 


In  the  class  ajfl*  >  we  must  draw  a  tangent  line  from  a  to 
b2  initially.  The  point  of  tangency  on  Bj  can  be  found  by 
using  the  formulas  developed  in  Chapter  V  of  this  paper.  We 
then  travel  around  the  obstacle  until  a  line  tangent  to  Bj  can 
be  drawn  to  z  (Figure  B.8). 
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a 


"2'  ^  al 

Figure  B.8   The  Shortest  Path  in  the  Class  a,^ 


Throughout  this  process  the  lengths  of  the  line  segments 
and  the  arc  lengths  representing  the  distance  traveled  around 
the  obstacles  is  summed.  The  final  sum  is  the  length  of  the 
true  shortest  path  for  the  respective  class. 

If  the  shortest  path  in  a  class  is  smaller  than  the  lower 
bound  on  the  next  class  to  be  considered,  then  the  search 
stops.  In  this  example  the  search  stops  after  the  first  class 
is  checked  since  the  length  of  aia^  ^s  1-8793  and  the  length 
of  afii    is  1.9032. 
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